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Abstract. In this review we consider those processes in condensed matter 
that involve the irreversible flow of energy between electrons and nuclei that 
follows from a system being taken out of equilibrium. We survey some of 
the more important experimental phenomena associated with these processes, 
followed by a number of theoretical techniques for studying them. The techniques 
considered are those that can be applied to systems containing many non- 
equivalent atoms. They include both perturbative approaches (Fermi's Golden 
Rule, and non-equilibrium Green's functions) and molecular dynamics based (the 
Ehrenfest approximation, surface hopping, semi-classical gaussian wavefunction 
methods and correlated electron-ion dynamics). These methods are described 
and characterised, with indications of their relative merits. 



1. Introduction 

This review is about the transfer of energy between electrons and nuclei. We focus 
especially on the theories and computer models that can be used to describe this 
process. Energy transfer is such a general concept, however, that we need to specify 
what exactly we have in mind here. A possible sequence of events that illustrates the 
kind of processes we have in mind is: 

(i) A system begins in equilibrium. 

(ii) It is then subjected to some external influence. This might be an electromagnetic 
pulse that transfers energy to the electrons thousands of times more efficiently 
than to the nuclei, thus giving the electrons a relatively higher energy. Or, at the 
other extreme, it might be a flux of neutrons which interact only with the nuclei 
which would have the effect of raising the energy of the nuclei relative to the 
electrons. (It could also be a number of other things, such as a beam of electrons 
or a flux of heat.) 

(iii) As a result of this departure from equilibrium, there will be a net flow of energy 
either from the electrons to the nuclei, or from the nuclei to the electrons, as the 
system moves back towards equilibrium. 

The flnal transfer of energy is made possible by interactions that couple the electrons 
to the nuclei. Of course these same interactions lead to other forms of correlated 
motion between electrons and nuclei, such as the modification of the electronic band 
structure, the effective coupling between electrons that produces superconductivity, 
and the coherent transport of electrons by virtual polarons. Interesting though they 
are they are not covered here. 

This review is structured as follows. In the next section we briefiy summarise some 
relevant phenomenology. We then describe the most important methods currently 
available for modelling these phenomena. Finally we briefiy consider the future of the 
theories. 

2. Phenomenology 

In this section we survey some phenomena observed in solids that are produced by the 
irreversible fiow of energy between electrons and nuclei. In what follows, neither the 
range of phenomena nor the list of citations are exhaustive. Instead we have merely 
attempted to provide a basis for thought and analysis for the interested reader. 
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2.1. Joule heating 

Probably the most familiar and easiest to observe phenomenon involving the inelastic 
transfer of energy between electrons and nuclei is the increase in temperature of an 
electrical conductor when a current passes through it. The idea that electric currents 
do work on conductors, and so cause heating, has been known since the mid nineteenth 
century when it was investigated by James Joule [1]. The treatment of transport 
in terms of Boltzmann's equation has been brought to a high level [2], with the 
coupling between the electrons and phonons treated perturbatively. This is accurate 
for macroscopic materials because no coupling to any one vibrational mode is large. 

With the rise in interest in mesoscale and nanoscale systems additional quantum 
effects became apparent, introduced by the small length scales. One such system 
is the atomic scale contact, which can be produced either by a break junction or 
a scanning tunnelling microscope. When a current is passed through the junction 
heating takes place that increases with applied voltage. Experimental evidence for this 
heating comes from observations of the voltage dependence of two level fluctuations 
and hysteresis at conductance discontinuities in atomic scale contacts [3,4] and from 
measurements on current induced rupture of atomic chains [5] . This heating, which can 
be very substantial, may at first seem mysterious since the electron mean free path is so 
much greater than an atomic spacing. However this net heating is simply a compromise 
between the small probability for an individual electron to scatter off a phonon in an 
atomic wire, and the huge density of current carrying electrons that accompanies the 
current densities attainable in quasi-ballistic metallic nanoconductors [6-8]. This can 
result in the apparent paradox of a hot conductor that is still largely ballistic, and 
thus resistance free, as far as individual electrons coming through are concerned. 

2.2. Relaxation of excited electrons 

The transfer of energy is not the only important consequence of the non-adiabatic 
interactions between electrons and nuclei. For example, non-radiative processes in 
semiconductors contribute to trapping of free carriers which has an impact on the 
conductivity and optical properties [9,10]. Because the binding energy of deep centres 
can be much larger than typical phonon energies, the phonons can participate in 
the optical transitions (which provide the most direct information about the defect 
electronic properties). Thus non-radiative processes can be observed in optical 
absorption spectra of lattice defects [11]. 

Paramagnetic ions in insulators can undergo spin reversal [12] (the ions switch 
between Zeeman levels). Initially an electromagnetic field is applied. This takes 
the population away from equilibrium. The system can then return to equilibrium 
through interaction with thermal phonons [13], a process that can be observed by 
spin resonance experiments [14]. There is intense interest at present in similar 
phenomena in nanoscale systems such as quantum dots, as the quantum confinement 
produces discrete electronic states which modify the relaxation rates of electrons and 
holes [15-17] 

Excitons can be thought of as a bound pair of particles: a negatively charged 
electron and a positively charged hole. Excitons can further interact with each 
other and form pairs [18], or bind to defects [19,20]. However, since the hole is 
simply the absence of an electron, with the positive charge originating with the 
nuclei, it is possible for the two quasi-particles to combine and annihilate with the 
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release of energy. This energy can escape as a photon or one or more phonons. In 
semiconductors a phonon may be essential to the recombination process because of 
crystal momentum conservation, though the presence of defects that trap excitons can 
remove this constraint [19]. The angular momentum state of the exciton can frustrate 
its radiative decay [20] , which is important for the efHciency of organic light-emitting 
diodes. A simple counting argument suggests roughly a quarter of excitons are in a 
singlet state (and can undergo radiative decay) while the remaining three quarters 
are in a triplet state (and cannot undergo radiative decay). However, there appear to 
be significant exceptions to this rule [21]. Excitons in quantum wells have increased 
binding energy because of their confinement, though their interaction with phonons 
appears unaffected [22]. In multiple quantum well structures, excitons can become 
trapped at a number of sites with different energies. Phonons assist in the migration 
at low temperatures (making up the difference in energies between binding sites) , and 
at higher temperatures (>10K) produce thermally activated migration [23]. 

2.3. Inelastic electron tunnelling spectroscopy 

Inelastic electron tunnelling spectroscopy exploits the transfer of energy between 
electrons and nuclei to probe the vibrational spectrum of molecules by passing an 
electric current through them. A schematic of the experimental arrangement is shown 
in figure nj There are two metallic plates (or contacts) across which a bias is applied 
{eV = fiL — Mi?) where V is the applied voltage and fiL and fiR are the chemical 
potentials of the left and right plates respectively, as shown in figure . Between 
the plates is an insulating region (often an oxide) that contains the molecules whose 
vibrational frequencies we wish to probe. In general there will be an elastic tunnelling 
current (an electron arrives at the right plate with the same energy that it had when it 
departed the left plate) whose magnitude increases linearly with the applied voltage for 
low voltages. At low temperatures, and for sufficiently large applied voltage {eV > hujQ 
where ojq is the angular frequency for the vibrational mode with lowest frequency that 
can exchange energy with an electron) , a conduction channel will open up in which an 
electron can arrive at the right hand plate and enter an empty state with a diminished 
energy, the energy lost being equal to one quantum of vibrational energy for the 
molecules in the sample. By increasing the voltage, additional channels open up 
corresponding to higher vibrational frequencies, or possibly multi-phonon processes. 
At higher temperatures it is possible for the electron to gain energy from the molecular 
oscillations. 

Lambe and Jaklevic [24] present a theory based on Fermi's Golden Rule (see 
section 13.2. 2|l suitable for oxides in which the electron exchanges energy remotely 
from the oscillators through an electrostatic dipole interaction. They also report 
experimental results for a system consisting of Al and Pb plates sandwiching AI2O3 
doped with small organic molecules. The experiments were performed from below 
IK up to 300K, and inelastic resonances corresponding to molecular vibrational 
frequencies were detected as peaks in the second derivative of the current-voltage 
curve {d^I/AV^). Klein et al [25] used a similar technique, and found that different 
ranges of voltage probed not only the frequencies of the impurities but also of the 
contacts and the oxide. The technique for absorbing molecules was improved by 
Simonsen et al [26] who used the liquid phase of the molecules. This made it possible 
for them to investigate the vibrational frequency spectrum of much larger molecules 
including proteins and nucleotides. Kirtley and co-workers have made a series of 
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Figure 1. The arrangement for inelastic electron tunneling spectroscopy. There 
are two metallic plates, between which is an insulating layer. A voltage is applied, 
so the electron chemical potentials of the two plates become offset, allowing 
electrons to flow from filled states on the left to empty states on the right. 
In the insulating sample there are molecules that can exchange energy with a 
tunneling electron: the electron can excite a vibration of angular frequency uj in 
the molecule, and so lose energy tko. 



contributions [27-31]. They have investigated the effect of the choices of metal and 
oxide for the tunneUing device, and found that these can influence the measured 
spectra [27,29]. They have developed a theory of the intensities of the vibrational 
modes using a transfer-Hamiltonian formalism with WKB wavefunctions [28] . In this 
model the interaction between the tunnelling electron and an oscillating molecule is 
through the Coulomb interaction with the partial charges on the molecule. It has been 
appHed to CHaSO;^ chemisorbed on alumina [31]. Multiple scattering Xa calculations 
showed that this long ranged potential indeed dominates the measured spectra [30]. 

Persson and Baratoff [32] showed that the excitation and immediate de-excitation 
of a molecular vibration by a tunneling electron can give a decrease in the resonant 
conductance. A closely related problem is the inelastic current voltage spectroscopy 
of a ballistic atomic metaUic chain [33,34], or a single resonant molecule between two 
electrodes [35]. The opening of new inelastic scattering channels, as the bias matches 
the energies of various phonon modes in the wire, also leads to a suppression of the 
electronic conductance, since, starting from ideal transmission, the conductance can 
only ever go down, and the inelastic I-V features take the form of dips in d^I/dV^. 

Scanning tunnelling microscopes (STMs) have been used for inelastic spec- 
troscopy. Inspired by this, Gregory [36] created a very small junction by using two 
crossing gold wires with argon and carbon monoxide in between. Resonant peaks 
were observed from hydrocarbon contaminants, as was Coulomb blockade behaviour. 
Ho and co-workers [37-40] have used an STM to probe the orientation, motion and 
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vibrational spectra of small molecules on metal surfaces. Yu et al. studied a molec- 
ular transistor and found shifts in the resonances when the gate voltage was varied. 
With the experimental advances has gone developments in the theory. Lorente and 
Persson [41] used density functional theory and a perturbative implementation of 
non-equilibrium Green's functions (NEGF) to study C2H2 and C2HD on Cu(lOO), 
and Galperin et al. [42] used NEGF and a model Hamiltonian to study Hne shapes 
and Hne widths of lETS features. 

2.4- Friction force on high energy particles 

The bombardment of materials by high energy heavy particles occurs in a number 
of situations. Ions with energies of order 10 keV are used for ion implantation 
in semiconductors [43], while in hydrogen fusion powerplants materials of quite 
extraordinary resilience are required to surround the region containing the plasma 
because they are subjected to bombardment by highly energetic neutrons (14.1 MeV) 
and other particles as well as exposure to very high temperatures [44]. 

Ion implantation of semiconductors is well understood [11,43]. The type of 
scattering experienced by the implanted ion depends on the mass of the ion and 
the relative angle of its trajectory to the crystallographic axis. The distribution of 
implanted ions is a gaussian centred about the projected range. The projected range 
in turn is determined by the stopping power which has contributions both from the 
nuclei in the crystal and from the electrons. The stopping power (S) is defined by the 
the projected range (Rp) through Rp = J^'dE/SiE), and is a sum of the nuclear 
and electronic components {S = Sn + Se)- The electronic component has a very 
simple form, namely Se{E) = keVs, where ke is a constant (in silicon it has the value 
lO'^VeV /cm). As the stopping power has the dimensions of a force, and is proportional 
to the velocity of the incoming ion {v = y^2E/M), it is behaving as a kind of friction. 

There has been considerable interest in the properties of candidate materials for 
fusion powerplants under extreme conditions, and extensive modelling work has been 
performed [45] . Clearly the collisions between atoms to form cascades is an extremely 
important process, but the dynamics of the atoms is significantly modified by the 
response of the electrons to this violent motion. The dynamics of radiation damage 
is typically divided into three phases: the displacement phase, the relaxation phase 
and the cooHng phase. The details of the interaction between the nuclei and electrons 
need not be the same in each phase because of the large variation in the energy of 
both the nuclei and the electrons between the phases. Further, the manner in which 
electrons transport heat is also modified when they become highly excited [46], or the 
lattice highly disordered [47]. There is still some controversy about precisely what 
happens, but there are some general ideas which are clear. Highly energetic nuclei can 
give up substantial amounts of energy to the electrons [48] , and in so doing experience 
an effective friction force [49]. As the coefficient of friction depends on materials 
parameters such as electron-phonon coupling strength and electron heat capacity, it 
leads to different rates of cooling in different materials, which in turn can modify both 
defect formation and healing of damage [49]. 

There is a rather different type of process in which electrons (light particles) are 
used to bombard a metal to produce damage (displacement of heavy particles) though 
at quite a low level (2 x 10~^ to 2 x 10~^ displacements per atom). Subsequent 
measurement of the variation of resistance as a function of increasing temperature 
during annealing of the damaged material provides direct information about the 
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diffusion barrier heights of point defects, which would be difficult to obtain any other 
way [50]. 

3. Theoretical methods and their applications 

As all theories involve some kind of simplification of reality, they must all throw away 
features that are considered unimportant. However, there will always be cases where 
what is normally unimportant suddenly takes centre stage. When thinking about the 
interaction of atoms with each other, it is usually a very good approximation to assume 
that the electrons are so light, and thus move so fast relative to the nuclei, that over 
the time it takes for the nuclei to undergo a significant displacement (say lO^^^s), the 
electrons will have undergone so many collisions, both with nuclei (possibly as much as 
1000 times in a simple classical view) and with each other, that they will have settled 
down into a well-defined state (normally taken to be their ground state, or more 
generally that state which minimizes the chosen free energy). Thus, we can treat the 
state of the electrons as known once the positions of the nuclei are given. This is the 
Born-Oppenheimer approximation, whose chief characteristic is that it allows us to 
treat the electrons and ions separately, and is ubiquitous because it greatly simplifies 
the process of understanding matter at the atomic scale. 

However, the Born-Oppenheimer separation is an approximation, and there are 
phenomena it cannot describe, notably those that are the subject of this review. 
To understand the nature of this approximation, and why it is unable to describe 
the irreversible transfer of energy between electrons and nuclei, we need only note 
the following. Even if the nuclei start at rest, collisions with electrons will result 
in small energy transfers to nuclei (of order 0.1% or less of the electron energy per 
collision, classically). Once they have started to fiuctuate, the nuclei can transfer 
kinetic energy back to the electrons in individual electron-nuclear collisions. But this 
two-way inelastic energy transfer depends on the velocities of the colliding particles. 
This dependence lies beyond the Born-Oppenheimer approximation. 

Fortunately, the amount of energy transferred in one collision is very small. Thus 
the departure from the Born-Oppenheimer approximation is a rather weak effect and 
so it makes sense to use perturbation theory. The perturbation can be characterised 
by those contributions to the forces on the nuclei that involve two electronic states, 
corresponding to electrons scattering from one state to another. The unperturbed 
electronic states are some representation of the Born-Oppenheimer energy surfaces, 
and for the nuclei harmonic oscillator states associated with small displacements from 
equilibrium positions on those surfaces are usually used (see [Appendix A[ . Lowest 
order perturbation theory can then be used to evaluate the rate of change of some 
quantity (such as the degree of excitation of the nuclear vibrations) as a result of 
the perturbation. For systems in which well-defined reference states can be defined, 
this prescription provides not only insight into mechanisms, but also numbers that 
can be compared with experiment. Examples include Joule heating in metals and 
non-radiative transitions at point defects in semiconductors. 

However, there are systems in which it is difficult to set up the reference systems. 
If the nuclei undergo substantial displacements, so that the concept of a reference 
position is not well-defined (for example, fiexible polymer strands in contact with a 
heat bath, or a metal subject to high-energy bombardment) then neither the electronic 
nor nuclear states are straightforward to characterise. However, such systems have 
the feature that much of their behaviour can be well described by molecular dynamics 



The transfer of energy between electrons and ions in solids 



8 



carried out within the Born-Oppenheimer approximation. So, the natural consequence 
is to see if the standard molecular dynamics algorithms can be revised to include the 
effects of the breakdown of the Born-Oppenheimer approximation. Since the full 
coupled problem is insoluble except in the simplest of cases we have to investigate 
possible approximations. 

In the following sections we discuss possible mathematical methods, both 
perturbative and molecular dynamics based, that can be used to investigate non- 
adiabatic processes. We begin with a simple classical model in which both the nuclei 
and the electrons are treated as classical particles obeying Newton's equations of 
motion. This cannot be considered a particularly realistic model, but it is simple 
to understand and in fact reveals much of the physics. We then proceed with the 
methods that make expHcit use of Born-Oppenheimer surfaces and oscillator states, 
both perturbative and non-perturbative. Finally we survey the methods based on 
molecular dynamics. 



3.1. Classical model 

While a classical model of the motion of electrons may not provide us with a 
quantitative scheme in general, it makes it easy to understand inelastic processes. 
To illustrate the points made we will look at the heating of an atom by a current of 
electrons [51] . 

The physical setup is as follows. We have an atom in a solid for which the 
influence of the neighbouring atoms is represented by a spring (an Einstein model) . 
The electrons we treat as independent particles which can interact with the oscillating 
atom. The Hamiltonian for the system is the sum of an electron Hamiltonian (He), 
an atom Hamiltonian (Ha) and a coupling Hamiltonian {He a) 

^ - E (ll + ^ir.)) + + \kx^) - E ^ • Mn) (1) 

i ^ ^ v_; J i 

where v{fi — X) is the interaction potential between electron i and the atom, and we 
have made the approximation v{ri — X) w v{ri) — J2i X ■ Vv{ri), which corresponds 
to small atomic displacements. Using Hamilton's equations we can write down the 
equations of motion for the electrons and for the atom 

^ or^^ 

mr,^ = ^ + y X,, ^'"^^''^ (2) 

In the absence of the coupling Hamiltonian IleA , from Eq. |2l we can write down an 
explicit solution for the atomic motion {X^{t) = sm{ ^/K/M t + (pi,)) . Furthermore, 
the electrons move independently of one another, with each electron moving with 
a constant energy (the scattering from the atom is elastic). Once we include the 
coupling Hamiltonian, this all changes. We can no longer write down a simple closed 
form expression for the atomic motion as it now depends on all the electrons; the 
electrons are no longer completely independent because their motion depends on X, 
which in turn depends on all the electrons; and electrons no longer move with constant 
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energy because they experience an external time varying force. It is the third point 
that is the main topic of this review. 

We can solve these equations for the case of motion in one dimension in which 
the atom appears as a hard wall potential to the electrons. In this case, by conserving 
energy and momentum during a collision, we get the following approximate expression 
[51] for the power delivered to the atom (w) by a current (j) of electrons in the Hmit 
that the electrons are much lighter than the atom 

«'«4j- (^) {K,~2Ka) (3) 

where i^eand Ka are the mean kinetic energies for the electrons and atom respectively. 
We see that there are two contributions, one positive from the electrons (which 
provides heating) and one negative from the atom (which leads to cooHng). The 
energy scale relevant to the electron kinetic energy is eV, where V is the appHed bias. 
At low temperatures, the energy scale relevant to the atom kinetic energy is hw where 
uj is the vibrational frequency of the atom. Thus for very small bias we might expect 
cooling, but for typical voltages we expect the electrons to heat the atom, as is indeed 
observed. 

Equations IS enable us also to see expHcitly the microscopic correlated electron- 
nuclear fluctuations that nucleate non-radiative inelastic processes and mediate the 
inelastic exchange of energy between the two subsystems. Consider an electron bound 
in an orbit with radius r and speed v around a nucleus. The nucleus starts off at rest, 
at the bottom of the parabolic well provided by the rest of the lattice. Let us solve 
Eq. 121 approximately as follows. To lowest order, ignore the motion of X in the second 
equation. Then the electron remains in the original orbit and provides a sinusoidally 
varying external force on the nucleus, of amplitude F = mv^ /r and frequency uj ~ v/r, 
along each of the two axes in whose plane the orbit lies. The nucleus in turn is now a 
driven harmonic oscillator. If we solve the equation motion for the nucleus assuming 
that it starts at rest, then the average kinetic energy over one period of oscillation of 
the nucleus is given by 



2 M 2 (tj2 _ t^2)2 

If a; ^ Wo, we see a resonance, with violent heating of the nucleus. This resonance 
is the classical analogue of the quantum transitions that become activated when the 
two frequencies match. In the Hmit lu^ u^, on the other hand, Ka settles at 

Ka = Kp 

where = mv^ 12. The nucleus has acquired some kinetic energy, and thus is no 
longer sitting still. The underlying motion of the nucleus, responsible for Ka above, is 
the wobble of a heavy particle induced by the passage of a light particle, famihar from 
planetary physics. As a consequence of this correlated wobble, nuclei are never truly 
frozen: momentum conservation, with the consequence that light and heavy particles 
always orbit a common centre of mass, leads to a repartitioning of energy between the 
two systems, in a ratio controlled by m/M . 

3.2. Born-Oppenheimer based electron-phonon coupling formalism 

In this section we consider those methods that make explicit use of energy surfaces, 
and the non-adiabatic coupling between them. 
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Figure 2. The parameters characerising a two level system. The Landau-Zener 
theory concentrates on the crossing region on the right-hand side of the diagram. 

3.2.1. Landau-Zener theory There is a generic model that has been successfully 
appHed to a huge number of problems: it is the two level harmonic model (or 
spin-boson Hamiltonian) and is often considered in the presence of a dissipative 
bath [52]. It has the advantage of providing useful insights pictorially as well as 
through numbers. In figure [2 is shown a representative system. The two diabatic 
surfaces (see [Appendix A| ) have the same curvature (vibrational frequency), but their 
minima are offset relative to each other both in energy and position. There are 
three other characteristic energies, two associated with vertical (for example, optical) 
transitions {Ei + Shw and E2 + Shu, where S is the Huang- Rhys factor), and one 
associated with non-radiative transitions {E2 + Eb). It is the latter that interests us 
here. 

The quantity Eb is the barrier height that a classical oscillator on the upper 
energy surface has to overcome to cross onto the lower energy surface. The transition 
between the surfaces is a quantum process (involving an electronic transition). This 
can correspond to a number of physical phenomena, with an important one in the solid 
state being electron capture (or emission) by a charge trap in a semiconductor [53], 
and another being atomic collisions with crystal surfaces [54]. This combination of 
quantum and classical descriptions was turned into a transition rate in 1932 both by 
Zener and Landau [55,56], and has been extended in a number of ways since (see for 
example [57,58]). 

The essence of the original Landau-Zener theory is as follows. The oscillating 
atoms (represented by just one coordinate, X in figure|3 and from now on referred to 
as the oscillator) move on a diabatic energy surface. The coupHng to the neighbouring 
surface is so weak (because of charge rearrangment in the case Zener studied) that a 
transition is unHkely, so we can neglect the influence of the transition on the dynamics 
of the oscillator for the purpose of determining the transition rate. No transition is 
possible except very near to the crossing. In this region the energy surfaces are treated 
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as linear, and the velocity of the oscillator as constant. By solving the time dependent 
Schrodinger equation for the electrons subject to the time varying potential due to 
the oscillator, the probability of a transition during one pass of the oscillator over the 
crossing is found to be [55] Px = cxp(— 2717) with 

^ 1^12 P f d{Ei-E2) 

^ h \ dt 

where E12 is the electronic matrix element coupling the two surfaces, and Ei — E2 
is the energy difference between the surfaces. The transition rate is then given by 
r = 2fPx [59] where / is the vibrational frequency. 




3.2.2. Perturbation theory There is a long history of using perturbation theory for 
studying non-adiabatic processes, both in the soHd state and elsewhere. A very 
important context for its use is in the Boltzmann equation treatment of the transport 
properties of solids [2]. Fermi's Golden Rule gives the transition rate from state \i) to 
state I/) 



r 



2 



{f\H,\2) 5{Ef-E,) (5) 



where Ei = (i|_ffo|*) and Ef — {f\Ho\f), Hq is the unperturbed Hamiltonian, and 
Hi is the perturbing Hamiltonian driving the transition. For the case of interest 
to us (electrons interacting with mobile nuclei) from [Appendix A| we can identify |z) 
with l^'nAf), I/) with l^'n'AT'), Ei with UnN and Ef with Un>N', where n indexes 
electronic states and N indexes nuclear states. To make the integrals in the energies 
and transition matrix elements tractable, the harmonic approximation is used for 
the nuclear states. For the adiabatic representation (see [Appendix A| | some further 
decisions about the dependence of the electronic wave functions on nuclear coordinate 
is also necessary. These decisions are implicit in the static lattice representation. 
In the context of non-radiative transitions the choice of approximation is discussed in 
detail by Stoneham in review articles [9,60] and a book [10]. Once the rates are known, 
experimentally observable quantities can be determined such as vibrational lifetimes 
of adsorbed molecules [61], or the power supplied to atoms {w) during electric current 
induced heating in nanocontacts [8,62]. The power is given by 

^ = ¥ II l(*n'W'|i?l|*nJV>|' {Wn' - WN)S{Un'N' " U^n) (6) 
n'N' 

where Wn is the energy of the nuclei in state N. These matrix elements have 
been calculated, within tight binding and density functional theory, for atomic and 
molecular wires [62-68]. The power into a phonon mode has the structure 

^ = S E /"(I - M l(/3|.9l«)l' {{N + ~ (ep + hco)] - NS[e^ - (ep - hco)]} (7) 

where |q;) is a single particle electronic orbital with eigenvalue and occupancy /„, 
and g is the gradient of the electronic Hamiltonian with respect to the oscillator normal 
coordinate. The phonon frequency is oj and its occupancy is N. For simplicity we have 
assumed that all comic masses are the same, M, but in practice this assumption is 
easily lifted. Equation accounts for the observed heating in nanocontacts [8,62,65]. 
The same perturbative model can also account for some of the inelastic current- voltage 
features in atomic and molecular wires [63,66-68]. 
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Lowest order perturbation theory clearly cannot handle multiple coherent 
transitions or strong interactions [69]. However, higher order theories can be used in 
these cases Csee 13. 2.311 . Problems where higher order theories have been used include 
the non-adiabatic transition at a crossing between two electronic levels in the presence 
of dissipation [58], cross-sections for inelastic scattering of electrons from a surface [70] 
and inelastic tunnehng cross-sections [71]. 



3.2.3. Non-equilibrium Green's functions Green's functions are one of the well 
estabHshed tools for deahng with problems in condensed matter. They provide 
information about the response at any point in space-time due to an excitation (usually 
the creation or annihilation of a particle) at any other point. Green's functions are used 
in perturbation theory, linear response theory, and quantum kinetic equations. Their 
real power becomes evident when we consider a system of many particles (electrons 
and nuclei) that interact, and more especially when the system is being driven out of 
equihbrium by some external forces. Their application to transport, in the framework 
of the Keldysh formalism, has a long history [72-78]. In the absence of particle 
interactions, the many-body Keldysh formalism reduces [79] to the steady-state 
transport formalism based on one-electron Lippmann-Schwinger scattering theory (for 
a brief review see [80]), which in turn reduces [78] to the Landauer formalism. In the 
following, we present some appHcations of the Green's functions technique for systems 
with interactions between electrons and atomic motion. Note that in this section (and 
the corresponding appendices) we depart from the usual notation for representing all 
operators by means of hats (such as H^) to avoid a confusion of extra symbols. Instead 
hats are used to represent only operators in the interaction picture. 

Beyond perturbation theory: Let us start by generalizing Fermi's golden rule to 
include higher-order terms in the perturbation V. The formal derivation based on 
scattering theory is rather involved, therefore we consider a simpler approach following 
the prescription given in Refs. [81-83]. This will permit us to introduce the different 
quantities that are needed in the full description based on non-equilibrium statistical 
physics [84,85]. 

Let us assume that the perturbation V in the full Hamiltonian H = Ho + V is 
turned on adiabatically in the distant past, i.e. V = Ve'^*'/^ (with 77 a small positive 
constant). Initially the system is, at time to, in the state \i) and evolves at time t > to 
into the state \i{t)): 

\i{t)) ^ c-'""'^''Uit, to) c'""'"/" \i) , (8) 

where U{t,to) is the time evolution operator (given within the interaction picture) of 
the system. The probability P{t) to find the system in the final state |/) at time t is 
given by |(/|i(i))P ^^d its time derivative gives the transition rate T fi = Ti-,f. To 
derive the generalized Fermi's golden rule, we include all orders of the perturbation V 
in the time evolution operator U. Then, we have 



with 



U{t, to) = 7^ I dti / ' di2 / ' di3 . . . / " ' dtS{ti)V{t2) . . . F(t„) c" 

„ U'^j Jto Jto Jto Jto 

°o 11/"* /■* / X //■/■* 
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where the perturbation V{t) is taken in the interaction representation V{t) ~ 
^iHot/hy ^-iHot/h ^j^j rp^ denotes the time ordering operator. Performing the time 
integrations and making use of the interaction picture for the perturbation, we obtain 

2 



pit)^\{fm) 



l2 



pTit/h 

-ifm 



(11) 



E., -Ef + iri 

and identifying the time derivative of P{t) with the transition rate T fi, one finds the 
generaUsed Fermi's golden rule: 

rf..^-\{f\m\'S{Ej^E,). (12) 

The quantity T in Eq. (Illjl is called the T-matrix and is given by the series expansion 

T^V + V V + V V V + ... (13) 

For the lowest-order expansion T ^ V ox equivalently by taking U(t,tQ) w 1 — 
i/h J^^dt'V{t') in Eq. lfT?Hl . one recovers the transition rate Tfi = Fi-*/ given by 
the usual Fermi's golden rule Eq. El 

The T-matrix can be rewritten in a compact and closed form 

T=V + V I T = V + VGIT (14) 
— Jio + 177 

which generates the same series as in Eq. El This allows us to introduce one of 
the different Green's functions, namely the retarded Green's function Gq defined, in 
the energy representation, by Gg(w) = [uu — + ir?]"^. This Green's function is 
obtained from the unperturbed Hamiltonian Hq. The retarded Green's function G'' 
for the total Hamiltonian H ~ Hq + V is defined by G^{uj) = [u — H + ir]]~^ . Then, 
the T-matrix can be expressed as T = -I- VGq T — V + VGV by making use of 
the Dyson equation that links the full Green's function to the unperturbed Green's 
function: G*" = Gq + GJjV^G'' = Gq + GqTGq. Now we have all the ingredients to 
develop the non-equilibrium theory, namely the concept of Green's functions, the time 
evolution operator and the Dyson equation. One already sees the importance of using 
Green's functions in determining transition rates when going beyond perturbation 
theory. Obviously, the conventional results of perturbation theory are obtained by 
expanding the Green's function in the Dyson equation or the T-matrix to the lowest- 
order in the perturbation. 

The T-matrix formalism with scattering boundary conditions has been used to 
study the effects of electron-atomic vibration coupling in the transport properties of 
atomic and molecular wires. This was done (i) to the lowest-order in the interaction, 
i.e. perturbation theory, in atomic wires [62-64] and molecular wires [65-67]; and 
(ii) to all orders in the interaction in model systems and molecular wires by using 
conventional Green's functions techniques [86-97]. The connection between the in- 
elastic T-matrix scattering formalism and non-equilibrium Green's functions has been 
discussed in Ref. [98]. 



Interaction as self-energy in the Green's functions: To describe the ground state 
or the thermal- averaged properties of a system of many particles, such as the coupled 
electron/atom system considered in this paper, it is useful to work with a many-body 
approach based on Green's functions. Such a Green's functions approach is even more 
useful and powerful when we want to calculate the properties of the system driven out 
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of equilibrium by some external forces. The exact definitions of the different Green's 
functions and their interrelation is given in [Appendix C| as well as the principles used 
to derive the theory for non-equilibrium conditions. 

Now, it is important to define precisely what the perturbation V means physically. 
Here we consider V as being a perturbation on the reference system (described by 
the unperturbed Hamiltonian _ffo) due to some kind of interaction. It is convenient 
to distinguish between different kinds of interaction, though generally speaking the 
calculation of such interactions and their inclusion in the Green's functions are formally 
equivalent for all kinds of interactions. In the following, it is convenient to consider 
as separate the following interactions: 

(i) The interactions between particles, i.e. the interactions between electrons, 
between phonons, or the interaction between electrons and phonons (or atomic 
motion). In practice, such interactions are incorporated in the corresponding 
electron or phonon Green's functions under the form of a proper self-energy via 
a Dyson-like equation for the Green's function. In the literature, it is usual to 
talk about electron Green's functions being dressed by the phonons, and phonon 
Green's functions being dressed by the electrons. The corresponding self-energies 
are usually difficult to calculate exactly for a many-particle system. They are 
obtained, up to some order in the interaction parameters, via the use of a many- 
body Feynman diagrammatic perturbation theory [99,100]. 

(ii) The interactions of the particles with an external field. As an example, one 
can consider an electromagnetic field exciting the electronic system by inducing 
electronic transitions, or an external applied bias that drives a nanojunction out 
of equilibrium by establishing a steady-state current fiow. The corresponding hot 
electrons would then transfer energy to the phonon degrees of freedom via some 
interaction of type (i) mentioned above. 

(iii) It is often interesting to partition a system into its different constituent parts. 
For example: a crystal partioned into a locaHzed region around a defect and the 
rest of the crystal; a surface with adsorbate partitioned into the bare surface 
and the adsorbate; a nano-junction partitioned into three parts (a central region 
whose transport properties are studied and the two electrodes to which the central 
region is connected). In this case, the interaction between two different parts (I) 
and (II) of the system (for example the part of the electronic Hamiltonian that 
couples regions (I) and (II)) appears under the form of a self-energy in the Green's 
functions expanded onto the subspace of one of these regions. Often these self- 
energies are also referred to as embedding potentials [101]. These self-energies are 
usually practical to calculate, especially for systems treated within a mean-field 
approach and when there is no crossing of the interactions of type (i) between 
regions (I) and (II). 

We now consider that the system has either reached a stationary-state regime or 
more simply is at equilibrium. Then the Green's functions depend only on the 
time difference t — t' , and their Fourier-transform G^{uj) depend on only one energy 
argument u. In the presence of interactions, the electron Green's function C'" 
obeys a Dyson equation (either at equilibrium or in a non-equilibrium regime): 
C-'^iuj) = Go'"(w) + GS'"(w)S'-''^(w)G'-''^(u;) where Gq'" is the Green's function of 
the electronic system in the absence of the interactions, and E'"'° is the self-energy 
arising from these interactions. A similar Dyson equation relates the phonon Green's 
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function D^'"^ in the presence of interactions to the undressed phonon Green's function 
-Dq " (see |Appendix C| |. For the non-equiUbrium state, the Green's functions G^'^ and 
D'^'^ obey another kind of quantum kinetic equation (see below or [Appendix C| |. 

The difficulty is to calculate exactly, or as accurately as possible, the different 
Green's functions by including all orders of the different kinds of interactions. It is 
also a challenge to calculate them numerically because the self-energies J]^=(''.«.<.>) 
arising from the interactions are actually functional of the different Green's functions 
themselves: S'^(uj) = T,^[{G^{uj)},{Gq{uj)}] where G"^ (Gg) denotes a Green's function 
in the presence (absence) of the interactions respectively. 

One way of solving the problem approximately is to expand the Dyson equation 
in a Born series of the Green's function Go in the absence of the interactions, 
G = Go + GoSGo + GoSGoSGoSGo + then choose lower-order Feynman diagrams 
in the interaction parameters for the proper self-energy S (the Born approximation), 
and finally substitute the Green's function Go in the expression of the self-energy by 
the full Green's function G. In this way, one introduces a self-consistent scheme since 
the Green's function G both determines and is determined by the proper self-energy 
S. This approximation is usually known as the self-consistent Born approximation 
(SCBA). Physically it means that the interaction is treated to all orders but that 
only a limited number of elementary excitations are generated, i.e. those given by 
the many-body diagrams chosen for the self-energy. Within the SCBA some processes 
involving crossed diagrams are omitted (if they were not already included in the self- 
energy). 

Now that we have described the principle of non-equilibrium Green's functions 
and the ways to include the interactions, we consider an example of their use in relation 
to electronic transport through a heterojunction in the presence of an electron-phonon 
interaction. Then we will briefiy describe how they can be used to derive quantum 
kinetic equations as a generalisation of the Boltzmann equation for transport. 

Electronic current in the presence of interaction: Let us consider a scattering 
central region (a quantum dot, a molecular or atomic wire) in which there are 
interactions between particles, and which is connected to two (left L and right R) leads. 
The latter are described by two non-interacting Fermi seas at their own equilibrium 
and thus characterised by two Fermi distributions /l and /«. The electronic current 
Jl flowing from the left lead into the central region is then expressed in terms of three 
Green's functions (G'''° and G*^) of the interacting central region. In the stationary- 
state regime, the current Jl is given by [72,73,75,78,102]: 

Jl = '-^J du;Ti{fL{oj)TLiu;)[G'-{uj)~G''{u;)]+rL{u;)G<{u;)} , (15) 

where the trace runs over the electron states of the central region and is related 
to the imaginary part of the retarded (advanced) self-energy E^^"^ arising from the 
coupling of the central region to the left lead (self-energy arising from an interaction 
of type (iii) — see above) . These self-energies and other self-energies arising from the 
other kinds of interactions should be included in the calculation of G^^°'\ 

An expression similar to Eq. ^1 can be obtained for the current Jr flowing from 
the right lead into the central region by interchanging the subscript L ^ R in Eq. ^| 
For a current conserving system, one has Jl = —Jb- The famous result of Meir and 
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Wingreen [78] is obtained from the symmetrised current J ~ \ {Jl — Jr) 

J - dc^TrK/LHTLH - fR{Lo)TR{Lu)) [G'-H - G'^H] + (FlH - T^M) G<M} (16) 

By using the relationship between the diflerent Green's functions: G^—C^ = G^ — 
G< {c.f. |Appendix~Cl l , Eq-QUcan be rewritten as Jl = e/h J dw Tr{S]^ (w) G> (t^;) - 
S^(aj) G<(tj)}, where the self-energies due to the coupHng of the central region to the 
left lead are cx I/lTl and oc — i(l — fL)^L- The physical interpretation of this 
equation is more transparent than Eq. ^1 The first term gives the current fiowing 
through the left contact from the left electrode towards the central region, because 
it includes G> which gives information about the empty non-equihbrium states in 
the central region and the self-energy E^ which is proportional to the distribution of 
occupied states in the left lead. The second term gives the current flowing through 
the left contact towards the electrode because it includes G^ which gives information 
about the occupied non-equilibrium states of the central region, and the self-energy 

is proportional to left lead empty states. 

Now, to get physical results for the current we need to calculate the different 
Green's functions involved in Eq. El These functions obey Dyson-like equations 
Qr,a _ Qr,a ^ G^'"^Y,rMQrM £qj. Qa,r another quautum kinetic equation for G^'^ , 

i.e. G<'> = (1 + G''E'~)Go + S-^G") + G''E<'>G". The first term in the equa- 
tion for G^^^ comes from the initial conditions and is seen as a boundary condi- 
tion. It is generally omitted for the stationary-state since it acts only in the tran- 
sient regime. The approach is now very useful if we can construct the self-energy 
functional that include the physics of the many-body system considered and that a 
sufficiently good solution can be obtained from the coupled equations for G"''' and 
G^'^ {c.f. [Appendix C| ). One important point is that the self-energy functionals 
and the solution obtained for the Green's functions preserve the condition of cur- 
rent conservation. It can be shown that within the decomposition made above for 
the different kinds of interaction, the condition for current conservation requires that 
/ dwTr{Er^gj(w)G>(w) - E.> gj(w)G< (w)} = 0, where Er^> are the self-energies due 
to interactions of types (i) and (ii) as defined above. This condition says that what 
is inelastically scattered into the central region should compensate what is inelasti- 
cally scattered out of the central region. It is not obvious that all choices for the 
self-energy functionals fullfil this condition. A good example is the SCBA for the 
electron-phonon interaction calculated from the undressed phonon Green's function 
Do- The SCBA has been used recently to study the effects of the coupling between 
electrons and atomic motion in atomic wires [103, 104]. Other self-energy functionals, 
developed in the spirit of the SCBA and including more elaborate approximations for 
the phonon propagators, have also been used to study the efltects of the interaction 
between electrons and atomic vibrations in model systems [42,105-110] and in more re- 
aHstic atomic and molecular wires [111-116]. However the problem of non-equihbrium 
transport through a coupled electron-phonon system is very complex to solve exactly, 
and more work towards this direction needs to be done. 

Quantum kinetic equations: It is possible to derive a quantum analog to the 
Boltzmann semiclassical equation for transport from the non-equilibrium Green's 
functions. The detailed derivation is rather complex and has been well described 
in other review articles [100,117,118]. Here we summarise the principles for obtaining 
such quantum kinetic equations. 
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We start with the lesser Green's function given in a space-time representation 
(r, t) for a fermion field operator in a solid. Then we perform a transformation 
similar to that given in [Appendix B| by defining the centre-of-mass and time-average 
coordinates {R,T) = i(ri + r2,ti + ^2) and the relative coordinates {x,9) = (ri — 
r2,ti — t2). The Green's function G^(r,0;R,T) expressed in these new coordinates 
is Fourier transformed in space and time with respect to the relative coordinates to 
give G^(fc,w;i?, T), from which is defined the so-called Wigner distribution function 
f{k,uj;R,T) — iG^{k,uj;R,T). The distribution f{k,u};R,T) is the quantum analog 
of the semiclassical distribution function f{r,v,t) used in the Boltzmann equation 
for transport (v being the particle velocity at point (r, In a system governed 

by the laws of quantum mechanics, because of the uncertainty principle and because 
scattering smears out the energy and momentum states, the usual relation between 
energy and velocity does not hold any more. Thus we have to consider a distribution 
function of the type f{k, u; R, T) where the energy to and momentum k are treated as 
independent variables. 

In order to get the quantum Boltzmann equation, we write the equation of 
motion for the Green's function G'^{k,uj; R,T). This is the most difficult part of 
the problem especially when we want to treat a many-body system with interactions 
[117, 118]. However the quantum distribution function /(fc, lu; R, T) is the only correct 
distribution function to describe particles in an interacting, many-body approach. In 
the presence of an electric field, the following quantum Boltzmann equation can be 
derived [117] 

where F is the force (electric field) acting on the electron. An additional driving 
term {v ■ Fdf/duj) is obtained in comparison to the usual Boltzmann equation. The 
last term {df/dt)s arises from the inelastic scattering due to interaction between 
particles. It can be calculated by the use of self-energies S(w) (arising from electron- 
electron or electron-phonon interactions) in the corresponding equation of motion for 
G^(fc, cj; i?, T). Once Eq. solved for the distribution function /(fc, u;; i?, T), we 

can calculate the macroscopic quantities, as measured in various experiments, such 
as the electronic density n{R,T), the electronic current j{R,T), the energy density 
nE{R,T) and energy current jsiR^T). The latter are obtained from the different 
moments (in energy uj or momentum k) of the distribution function /. Finally, the 
semiclassical Boltzmann distribution f{R, v,T)is found by integrating the distribution 
function /(fc = mv, lu; R, T) over energy. 

3.3. Molecular dynamics based methods 

There have been many suggested ways to modify the standard Born-Oppenheimer 
based molecular dynamics algorithms. It is usual to speak of their being broad 
categories of methods: surface hopping and effective interaction. The former category 
is very popular for studies of molecular systems, but has well-known deficiencies (it 
is somewhat ad hoc, and coherence between trajectories is lost). The latter approach 
(usually in the Ehrenfest approximation) is popular for some problems but also has 

I Note that the function /{fc, u); R, T) = iG< (fc, ui; R, T) is related to the Wigner matrix used in CEID 
(see sention l3.3.4^ by an integral over the energy u! [119]. However, as used in CEID the transform is 
over nuclear degrees of freedom, while here it is over electronic ones. 
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limitations. One recent method (Correlated Electron-Ion dynamics [120-123]) sits 
somewhere between the two usual categories, though it is viewed as an extension of 
the Ehrenfest approximation. We begin with a discussion of the Ehrenfest method 
because of its importance as a starting point for the other methods. We then present 
alternative methods that attempt to overcome its deficiencies. 



3.3.1. Ehrenfest method This is a simple, but ingenious, method that successfully 
adds some non-adiabatic features to molecular dynamics, but is incomplete. There 
are many ways to derive the formal equations for this method, but the one we shall 
use here most naturally leads us to CEID which is discussed in detail below. 

The Ehrenfest method consists of two separate approximations: the electrons 
interact with the nuclei through a classical field generated by its charge distribution 
(a mean field approximation, and the main source of error); the nuclei obey classical 
(Newtonian) mechanics. We will work with the density matrix. 

The first approximation allows us to write the full density matrix p as a tensor 
product of an electron density matrix pe and a nuclear density matrix p^: p = pe®PN- 
If we substitute this into the quantum Liouville equation 

and then trace out either the nuclear or electronic degrees of freedom we obtain the 
following equations of motion: 

'-^-kin.M (19) 

where = Tr^r |iJePAr| , = Tjs + Tre |//e/5e|, Tr^v {• ■ •} means a trace over 
nuclear coordinates, and Tre {• ■ •} means a trace over electronic coordinates. The total 
Hamiltonian {H) has been separated into the nuclear kinetic energy (T/v) and the rest 
{He)- It is the presence of the traces in the definitions of the effective Hamiltonians 
that results in the electrons and nuclei only responding to the classical fields. Note that 
below we show how the Ehrenfest approximation can lift the mean field approximation 
for the nuclei because they can be treated as classical particles with unique trajectories. 

The classical nature of the nuclei is introduced by replacing the quantum 
commutator brackets by classical Poisson brackets. The equation of motion for the 
classical nuclear density matrix (now a phase space density) becomes 

Apn{R, P) dHNjR) dpNjR, P) P. dpNjR, P) \ , . 

dt dR^ dP^ A'U dR^ J ^ ' 

where i?^ and P^, are the position and momentum coordinates and is a nuclear 
mass. If we consider just a single classical trajectory (indicated by the subscript T) 
then we have: 



PN 



{R, P) = Y[S{R,- RrAt)) 5 (Pi^ - PrAt)) (21) 



where RT,u{t) and PT,v{t) are the nuclear positions and momenta along the trajectory. 
Substituting Eq. [2] into Eq. 1201 we obtain 
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Inserting this nuclear density (Eq. |21 into the equation of motion for the electrons 
(Eq. we obtain 

f = rn[^Aip.] (23) 

Equations [22l and [231 constitute the Ehrenfest approximation. The use of a single 
nuclear trajectory in fact takes us halfway to introducing microscopic fluctuations: 
the electrons now respond to individual nuclei. However, the converse is not true, 
with the nuclei still responding to an average distribution of electrons. 

This approximation can fail either when the nuclei have to be treated as quantum 
particles (for example, when tunnelling takes place), or the nuclei respond to the 
microscopic fluctuations in the electron charge density as well as the mean (for 
example. Joule heating). In this review we are principally concerned with the latter 
case. 

To make this more transparent, let us suppose that the electronic density matrix 
can be represented by a phase space distribution pe{r,p) (this can be formally justified 
by means of the Wigner transform as described in [Appendix B| |. In this case we have 
Hn{R) ~ J drdppe{r,p)He{R-,r,'p)^ and hence the force on the nuclei {F^) due to the 
electrons and the other nuclei is (Eq. [22|l F„ — — J drdppe{r,p)dHe{R; r,p)/dRu. If, 
as for the nuclei, we could isolate a single electronic trajectory rrit), then the force 
would just be Fi, = ~dHe{R; rT,PT)/dR^. The integral, then, gives the average force 
produced by a set of trajectories with each one contributing to the force with a weight 
Afdppf,{r,p). 

This averaging has the effect of reducing the ability of electrons to pass energy to 
the nuclei. To understand this, consider the following simple model. Suppose we have 
one nucleus that experiences forces from the electrons that have two contributions: 
F = —kX + f{t). The first, harmonic, term corresponds to motion on a Born- 
Oppenheimer surface. The small residual non-adiabatic corrections are represented 
by the force f{t). Let us define the energy of the nucleus as Un = P'^ /2M + 

The rate change of this energy (the power given to the nucleus) satisfies = X ■f{t). 
If we now evaluate the average power supplied over some time interval r we obtain 



1 



to+T 



Un) =- X-f{t)At (24) 

' ^ ^ Jto 

If f{t) varies only slightly during a vibrational period of the nucleus, then the power 
will be very small, as the average velocity of the nucleus is zero. If it fiuctuates rapidly 

but there are no correlations between X and /, then the fiuctuations average out and 

the power again vanishes. However if we allow X to be a functional of /, via the 
equation of motion, then in Eq. there are correlated force- velocity fiuctuations, 
which can in turn be expressed in terms of the force-force correlation function [51]. 
It is these correlated fiuctuations that generate the power. This result is known 
as the fiuctuation-dissipation theorem. The effect of the averaging over electronic 
trajectories in the Ehrenfest approximation is to break the microscopic correlations 
between the force experienced by the nucleus due to the electrons and the momentum 
of the nucleus, and thus to suppress the energy that can be transferred to the nucleus. 
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This breaking of correlation is caused by the averaging of the electronic trajectories 
for every nuclear configuration. 

The above analysis is most natural for metallic systems which we can think of as 
nuclei embedded in a gas of rapidly moving light electrons. For small molecules with 
well spaced energy levels the criticism of the Ehrenfest approximation can be made 
differently. The calculations usually involve an electronic transition between two levels, 
induced by the momentum of the nuclei. It is found that the results are often in poor 
agreement with accurate quantum calculations. This is explained by noting that the 
single nuclear trajectory experiences forces simultaneously from a number of partially 
occupied energy surfaces, whereas it should be experiencing it from one only, but 
with separate trajectories on each surface. This is just another manifestation of the 
averaging discussed above. 

These problems manifest themselves in the failure to produce certain outcomes 
in simulations, such as heating of nuclei by current carrying electrons [51], thermal 
equlibrium between electrons and nuclei [124,125], and correct transition probablilities 
[126]. Various schemes have been proposed to overcome these problems. One starting 
point is to recognise that the ionic wavefunction rapidly decoheres after a transition 
leading to modified equations of motion [127-129], or to a stochastic algorithm in which 
trajectories on separate Born-Oppenheimer energy surfaces are treated as independent 
(see Section K^.^i.2l below). Another is to treat the electrons and nuclei as having 
correlated motion leading to wavepacket-like methods [130, 131] (see Section 13.3.41 
below). 

While the Ehrenfest approximation is often presented in the literature as a straw 
man (that is, a method used simply to show how much better other methods are), it 
is also used to obtain useful results. Two applications are: 

(i) Time-dependent density functional theory has been used successfully to study 
the friction forces experienced by hydrogen atoms scattering off a metallic 
surface [132]. The friction is due to the excitation of electrons by the moving 
atoms. Energy transfer in this direction is correctly described by the Ehrenfest 
approximation. 

(ii) Polymers have been successfully studied quite extensively using the Ehrenfest 
approximation within a simple tight binding model for the electrons. 
Investigations have included: polaron drift in an applied electric field [133]; the 
dynamics of polarons formed after the excitation of electrons by photons [134]; 
the nucleation of stable self-localised excitations [135]. 

3.3.2. Surface hopping methods Surface hopping was originally proposed by Tully 
and Preston [136] as a phenomenological extension of the classical trajectory method 
which incorporates non-adiabatic effects. In the classical trajectory method, an 
ensemble of trajectories sampling a series of initial conditions for a gas phase chemical 
reaction is integrated in time. Quantities such as reaction cross sections are obtained 
as averages over the ensemble. 

As mentioned in the previous section, within the Ehrenfest approximation nuclei 
feel a force from the electrons that is an average over many electronic trajectories. In 
cases where adiabatic electronic states are close to each other in energy over relatively 
small regions of space and are well separated otherwise, as in the case of a non- 
adiabatic chemical reaction, when the nuclei leave a region of strong coupling the 
force felt by them is a weighted average of the force in each surface. Tully and 
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Preston observe that "intuition requires that the motion take place on one surface or 
the other". That is, after going through a region where there is significant coupling 
between surfaces the ions must stay on the original surface on which they were moving 
or jump to a different one. This intuition is guided by the fact that in a molecular 
beam experiment, reactants convert to one set of products or the other, but one 
does not detect states which are mixture of different product channels. This implies 
some decoherence of any mixed state before the outcome is measured (or during the 
measurement) , which in turn suggests a picture of the dynamics in which nuclei move 
adiabatically unless they encounter a region of strong non-adiabatic coupling. This 
can be either a point where adiabatic surfaces cross or where they are close in energy. 

This idea was implemented into the first realization of the surface hopping 
procedure by determining a fixed region or set of points in configuration space where 
switching between adiabatic states is likely. When the system passes through those 
regions during the dynamics the decision to hop to a different surface is taken according 
to a probability calculated from the Landau-Zener expression (see section 1^. 2. After 
a jump to a different potential energy surface has taken place, the nuclei change their 
potential energy by a finite amount, which is small if the jump occurs when the two 
surfaces are close to each other. In order to conserve energy the surface hopping 
procedure must be complemented with a prescription on how to change the velocities 
of the nuclei after a hop. In their original proposal, Tully and Preston rescaled the 
velocities in one particular direction, chosen from the topology of the potential energy 
surface and the calculated couplings. As with the standard classical trajectory method, 
properties of interest are calculated from the properties of an ensemble of integrated 
trajectories. 

Tully later extended the method by lifting the constraint of fixed hopping regions 
and named the new method Molecular Dynamics with Quantum Transitions (MDQT) 
[137, 138]. In this method, the jumps between surfaces can occur at any point during 
the trajectory. The decision whether a hop should occur and to which surface is 
taken according to a probability of a jump occurring. This is calculated from the 
time evolving state populations, which in turn is obtained from the simultaneous 
integration of the nuclear dynamics on a certain adiabatic potential energy surface 
and the equation of motion for the electronic wave function or density matrix. This 
starts by recasting the equations for Ehrenfest dynamics given before in terms of the 
adiabatic electronic states (see section [Appendix"^ . First, the following anzatz is 
made for the electronic wavefunction 

*(t) = ^c,W$,(i?(i)) 

j 

where $„(-R(t)) are the adiabatic eigenstates of the instantaneous electronic 
Hamiltonian, and are a function of time via the nuclear coordinates. Replacing this 
ansatz in the Schrodinger equation for the electrons we obtain the following equation 
of motion for the electronic wavefunction coefficients: 

ihcj = CjEj {R{t)) - ift ^ CkR ■ dkj (25) 

k 

where the non-adiabatic coupling vector dkj is defined by 

4 = / <i>^(^(t))v^$,(i?(t))di?. 
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The same result can be obtained via a variational procedure from a suitably 
constructed Lagrangian [139]. The second term in Eq. |2ll is the one responsible 
for transitions between states. This transitions will be more effective if the nuclear 
velocity points in the direction of the non-adiabatic coupling vector. The rate of 
change of the population on a given state is given by the equation of motion for a 
diagonal element of the electronic density matrix: 



Where we have taken into account the fact that the matrix of non-adiabatic coupling 
vectors is anti-hermitian. The fewest switches algorithm proposed by Tully was 
designed so that the choice of where to hop resembles the kinetic Monte Carlo 
integration of a master equation for the state occupations given by Eq. [SEI In a 
given time step the probability that a trajectory occupying state i hops to a different 
state j is proportional to bij . This probability is constrained in a way that minimizes 
the number of hops required to achieve consistency between the occupations and the 
number of trajectories in the ensemble that occupy a given state. Tully demonstrates 
that the algorithm ensures that at any time the number of trajectories in a given 
state is representative of the occupation of that state. However, the problem with this 
argument is that there is no unique set of occupations, but instead there is one for 
each trajectory. If the occupations of the states from different trajectories diverge fast 
relative to one another because the nuclei follow different paths (something that will 
occur faster the more hops there are) then consistency will be lost [140]. 

In order to conserve energy whenever a hop occurs the velocity component 
pointing in the direction of the non-adiabatic coupling vector is rescaled. If, from 
the algorithm, it is determined that a hop must occur to a potential energy surface for 
which there is not enough energy, the hop is aborted and the velocity in the direction 
of the non-adiabatic coupling vector is reversed [138]. These classically forbiden hops 
can become a problem if they occur too often since they contribute to breaking the 
consistency between the propagated occupations and the number of trajectories in a 
given state [140]. 

Comparison with fully quantum results for simple systems shows that the method 
produces qualitatively correct results [137, 141]. Unlike the Ehrenfest method it 
reproduces the correct Boltzman populations when the system is coupled to a 
thermostat [125]. The fact that for each particular trajectory a full quantum 
wavefunction is propagated means that the method accounts in some way for coherence 
within each trajectory and some effects due to interference can be reproduced [137]. 
However no interference exists between trajectories. Further, it is difficult to asses, 
because of the way in which the interaction between nuclear and electronic degrees of 
freedom is treated, whether the resulting quantum coherence effects are meaningful 
[142]. One particular problem is the fact that the method seems to work only when 
the electronic states are represented in an adiabatic basis [143]. This is probably due 
to the fact that the use of this basis leads to the least number of hops [144]. The 
rescaling of the velocity after a hop, which is another controversial feature of surface 
hopping, can be justified in a number of ways by semiclassical arguments [145, 146]. 
Although the replacement of the propagation of a single, fully coherent, wavefunction 
by an ensemble of trajectories can be made plausible by semiclassical arguments [143] 
some aspects of the method are clearly ad hoc making it difficult to improve in a 
systematic way, or to predict in which situations it can be appHed successfully. 




(26) 
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Extensions of the method exist for cases in which there are both discrete and 
continuum states [147]. They either use different schemes to choose when to hop [148] 
or replace sudden hops by a continuous switching [149]. Prezhdo and Rossky proposed 
a combination of the Ehrenfest and surface hopping procedures in which the force used 
to propagate the nuclear trajectory comes from the integration of Ehrenfest equations, 
not just one adiabatic state, but this wavefunction is reduced to a single adiabatic state 
when a hop occurs or when the validity of the Ehrenfest approximation is violated 
according to some pre-established criteria [150]. 

Tully's MDQT is one of the most used methods to deal with problems in 
which coupling between electrons and ions is fundamental. In particular, it has 
been extensively used to study non-adiabatic processes in liquids, given that the 
complexity of liquid structure prevents application of other techniques. One of the 
first applications of the method was to the problem of the solvated electron. Space 
and Coker explored the relaxation of an excess electron in dense fluid helium. MDQT 
provided information on the dependence of the electronic relaxation process on the 
initial electronic state [151]. Later, MDQT and some of its variants were applied to 
an excess electron in water [152-154]. Another classical problem to which the method 
has been applied is the simplest photochemical reaction in solution, photodissociation 
of a diatomic molecule. Coker's group pioneered the application of MDQT to this 
problem, studying the photodissociation of I2 in liquid and solid rare gases [155, 156]. 
Others later studied photodissociation of diatomic molecules embedded in rare gas 
clusters [157,158]. The case of the ICN molecule which can isomerize after dissociation 
has been studied in rare gas matrices [159], where some of the vibrational modes of 
the molecule were included within the quantum description, in bulk water [160] and 
at the liquid/vapour interface of water [161]. An important effect of the solvent in 
the case of photodissociation is that some processes which occur with high quantum 
yields in the gas phase are inhibited in the solvent. This effect is explained in 
terms of the caging of the reactants which are held together in the solvent and 
have the opportunity to recombine. Particularly dramatic is the case of azomethane 
(H3CNNCH3) which photodissociates in the gas phase but isomerizes around the NN 
double bond in the solvent. This problem was studied using surface hopping in the 
group of Persico [162-164], becoming one of the flrst non-adiabatic simulations of 
condensed phase photochemistry with a realistic solvent and a polyatomic solute. 
Since photodissociation is a classic problem in photochemistry, there are abundant 
experimental results for these systems which can be compared with the results obtained 
from the simulations. In general the agreement has been found to be good. 

Relaxation of a charge transfer excited state via intramolecular electron transfer in 
solution is another classic example of a non-adiabatic problem in solution chemistry. 
Lobaugh and Rossky studied the relaxation of excited betaine in acetonitrile. The 
electronic structure of the molecule was described using configuration interaction 
applied to a simple semiempirical model [165]. The authors found that besides the 
motion of the solvent, which is effective during the initial stages of the relaxation from 
the excited state, also some internal motions of the molecule played a fundamental 
role in the relaxation process. 

Photoisomerization around a double bond is a paradigmatic example of non- 
adiabatic chemistry, and important biological functions such as vision hinge on it. 
Surface hopping has been used to study the photoisomerization of butadiene [166], the 
chromophore of the photoactive yellow protein [167] and of retinal in bacteriorhodopsin 
[168]. Azobenzenes are one class of compounds that can be used for the generation of 
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materials which change their properties via illumination. Trans to cis isomerization 
of azobenzene is the key for this effect, and surface hopping has recently been used to 
solve some controversy on the detailed mechanism of the process [169]. 

Very few applications of surface hopping to processes in solids exist. Yokozawa 
and Miyamoto studied the breaking of Si-H bonds from H terminated vacancies in 
Si02 by hot electrons using a primitive combination of the surface hopping procedure 
and first principles molecular dynamics [170]. Bruening and Friedman used surface 
hopping and a tight binding Hamiltonian to describe charge transfer from a conducting 
polymer to Ceo molecules [171]. Bach and Gross treated with surface hopping the 
problem of charge transfer in molecule-surface scattering [144]. 

An important difficulty with the implementation of surface hopping procedures 
is the need for accurate adiabatic states and non-adiabatic coupHng vectors which 
must be obtained from some model description of the electronic structure or as a 
parametrization of ab initio results prior to the simulation. Recently, surface hopping 
has been implemented in ab initio molecular dynamics methods with on the fly 
calculation of the electronic structure. Doltsinis and Marx have implemented surface 
hopping within a restricted open-shell Kohn-Sham scheme [172] which allows the 
calculation of the low lying excited states and the corresponding coupling vectors 
between states [173]. The method was applied to cis-trans photoisomerization of 
formaldimine (H2CNH) [173], excited state proton transfer and internal conversion of 
o-hydroxybenzaldehyde [174], and the photostability of methylated DNA bases [175]. 
Prezhdo et al. have implemented a similar scheme where the adiabatic states used 
for the surface hopping propagation are different excited Kohn-Sham determinants 
and applied the method to study the nonradiative relaxation of the chromophore 
of the green fluorescent protein and electron injection from alizarin into titanium 
dioxide [176]. 

Recently, some new formulations of surface hopping techniques have appeared 
which are constructed from well defined approximations of the exact quantum 
dynamics [177-179]. These techniques have up to now been demonstrated to produce 
excellent results when compared with exact calculations for model systems, and some 
small problems [178]. Perhaps in the future these formally correct methods will become 
available for realistic condensed phase systems. 

3.3.3. Frozen Gaussian Approximation and related methods In this section, we will 
cover the Frozen Gaussian approximation (FGA) and various related semiclassical 
(SC) techniques including the Initial Value Representation (IVR) and the Herman- 
Kluk propagator. We will also consider methods which use the FGA for nuclear 
wavefunctions (one of these, surface hopping, which sometimes uses FGA for the nuclei 
is discussed in Section f3.3.2|l . There are various other reviews of these techniques 
which go into more detail: a thorough, though very formal, review of semiclassical 
methods [180]; an overview of SC-IVR methods [181]; a discussion of wavepacket 
(FGA) and surface hopping methods [182]; and a more general review covering time 
dependent methods for large systems [183]. 

The use of gaussian wavepackets as a basis for nuclear wavefunctions started 
in scattering calculations which relegated the electrons to the role of providing a 
potential; their use was pioneered by Heller [184-189]. From this work, he proposed 
that a simple way to broaden single classical trajectories would be to use a Gaussian 
function whose width was fixed and whose average position and momentum followed 
that trajectory [190]. For a gaussian centred on the phase-space point (r, p) with 
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width 7, we write: 



x|r,p)=( — ) cxp[-7(x-r)^ 



N/4 



ip 



c)/h] 



(27) 



The semiclassical approximation (due originally to Van Vleck [191] as an extension 
of the WKB method to time-dependent problems) provides a way to address dynamics. 
Consider the amplitude to change from state 1 to state 2 for a system with continuous 
position and momentum (r, p): 
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The semiclassical approximation changes the propagator to: 
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5(r2,ri) 



dt' {T{t') - V{t')) 



Here S'(r2,ri) is the classical action for the trajectory going from ri to r2 in time t\ 
this makes clear the connection to Feynman's path-integral formulation of quantum 
mechanics. The factor C contains a number of factors which there is not space to 
discuss here (see for instance Ref. [181,192]). The sum over roots arises from the 
stationary phase limit [193] which gives all classical paths for which 5S = 0. Solving 
this equation as written gives a non-linear boundary value problem, where all values 
of pi which yield Y2 must be found; this will lead, in general, to multiple roots. 

The initial value representation replaces the integration over final coordinates, 
dr2, with integration over initial momenta, dpi (bringing in a Jacobian with the 
change of variables) . The sum over roots also disappears because the initial conditions 
determine a unique classical trajectoty. We can write: 
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This approximation allows some quantum interference and tunnelling effects to be 
described, while involving only real, classical trajectories. 

The most commonly used combination of a Frozen Gaussian basis, 
semiclassical IVR is the Herman-Kluk propagator [194,195]. 

f d^r^d^pi 

/ in \\N ' (xkt.Pt)C(r»,Pi,^) 



and the 



X exp[iS'(rj,pj,t)/fi] (rj,p,|'0o), 

for an TV-dimensional problem. A trajectory starts at phase space point (ri,pi) and 
runs for time i to a phase space point (rt, pt). The form of the semiclassical prefactor 
required is important if the calculation is to be carried over long time periods [195]: 

1/2 



C(r„p,,t) 



1 / 9pt dvt ■f-'^^t * ^Pt 



(31) 



The SC wavepacket methods described thus far have the important features of 
time-reversal and unitarity [196]. However, in their original form, there are often 
rapid oscillations in the integrand which cause problems with convergence. These 
oscillations were shown to be significantly damped [197] by the merging of the cellular 
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dynamics method [198] with the Herman-Kluk propagator. Other solutions to this 
problem include time integration over short periods [199] yielding effective averaging 
over the short period. Still within the realm of reactive scattering, the limitations 
of the Herman-Kluk propagator have been explored in 2D and 4D modelling of H2 
scattering from Cu(OOl) [200]. The 2D simulation was not sufficiently accurate, 
while the full 4D simulation required inordinate amounts of computational effort. 
There has been appHcation to non-adiabatic electronic evolution, in particular to the 
canonical problem of electron solvation [201]. In this area, FG propagation coupled 
with perturbation theory has been compared to the surface-hopping techniques [202]. 
In regimes where tunnelling is important, surface hopping is less accurate, though the 
FG propagation is only accurate over short timescales, and can depend on the width 
chosen for the Gaussian. Some analysis of the accuracy of using classical mechanics 
for the propagation of the wavepackets in condensed many-body systems [203] showed 
that while the wavefunction can spread (losing the justification for the FGA), densities 
continue to be localised, due to interference effects. These ideas have been applied 
to the calculation of non-linear optical response functions, with reasonable accuracy 
[204]. These SC-IVR methods scale exponentially with the number of degrees of 
freedom [205] though this does depend on the number of final states contributing and 
the energy resolution required. 

In order to model non-adiabatic transitions, the FGA (and any SC-IVR) must 
be extended. It has been shown [192, 206] that it is possible to linearise the SC- 
IVR [207], splitting the entire system into a system and a bath; and expansion is then 
made in terms of the bath coordinates, using Wigner distributions. This technique 
has been applied to simple models (a Morse oscillator coupled to a single harmonic 
mode) [192], and the spin boson problem with different spectral weights [206,208], 
with good agreement to exact results. 

The multiple-spawning technique [209-211] uses frozen Gaussians as basis 
functions during a time-dependent simulation. An effective non-adiabatic coupling 
is monitored between two states / and /'; using a diabatic (see [Appendix A| | 
representation for the electronic wavefunctions, this can be written as : 



d(R) 



(imn 



VH{R)-Vrp{R) 



(32) 



When this coupling for state / exceeds a threshold, new nuclear wavefunctions are 
added (or spawned) on the new state /' at a constant rate (so that more wavefunctions 
are added the longer the system remains in the area of non-adiabatic coupHng). 
These new functions are then propagated, and can themselves spawn. This leads 
to a system whose nuclear degrees of freedom follow multiple trajectories on different 
surfaces simultaneously. It has been applied successfully, among other things, to 
a two dimensional non- reactive collision [211], and more recently to light-driven 
reactions [212]. 

The Zhu-Nakamura (ZN) theory [57] for non-adiabatic transitions (an extension 
of the Landau-Zener formalism that works with adiabatic states) has been combined 
with the use of frozen Gaussians for the nuclear wavefunctions [213-215] to extend 
both theories, both with the surface hopping formalism [213,214] (specifically using ZN 
theory to calculate the non-adiabatic transitions) and as a separate method [215]. This 
approach has the advantage of giving the rates from an analytic theory, while allowing 
nuclear wavefunction propagation, and works well in multidimensional problems (by 
comparison to more complicated and costly numerical calculations of the problems) . 
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Figure 3. The central idea in correlated electron ion dynamics is that the width 
of a nuclear wave packet is narrow on the length scale of the separation between 
nuclei. That is, W <^ d where d is a typical distance between nuclei. In this figure 
the mean density of the nuclei as a function of position R is given the symbol pjv. 
Correlations between the electrons and nuclei are maintained by allowing the 
electronic state to vary as a function of position of the nuclei within the wave 
packet . 



3.3.4- Correlated electron-ion dynamics The small amplitude moment expansion is 
a scheme to introduce the correlations between the electronic fluctuations and the 
nuclei that are missing in the Ehrenfest approximation [51]. When this expansion is 
combined with molecular dynamics, we refer to it as Correlated Electron-Ion Dynamics 
(CEID). The starting point is the Ehrenfest equations (which are, of course exact, and 
which give the name to the approximation): 

Ri, 



P. 



Trip 



dHe 



(33) 



where Ry = Tr and P^, = Tr|p^p|, and Ri, and Pi, are the position and 

momentum operators for the nuclei. In the Ehrenfest approximation we were able to 
work with a single trajectory: in that case R^ = Rtm and P^ = Ptm- However, since 
the Ehrenfest equations refers to quantum nuclei, the uncertainty principle makes it 
impossible to have individual trajectories. The best we can manage are narrow wave 
packets (see flgureEJ. In this case Ri, and Pn correspond to the trajectory of the mean 
of the wave packet. 

If the wave packets really are narrow (which they often are in condensed systems) , 
then we can expand He{R) about R in the fohowing way [120] (where R and R stand 



The transfer of energy between electrons and ions in solids 



28 



(34) 



for the set of coordinates {i?^} and respectively): 
where Ai?^ ^ — R^. Then the mean force (Eq. I33II satisfies 



where pe = TrAr {p}, /ti,^ = TrAr |/5A^^|, and il2.uiy' = TrAr ^p/\R^/\Riy>^ . The first 
term in Eq. ISHis the Ehrenfest approximation. The higher terms account for the fact 
that the wave packet has some finite width, and that the mean force must include 
an average over the paths included within the packet. However, at this point it is 
important to note a fundamental distinction. There are two separate contributions 
to the total width of the wave packet. The more obvious is just the normal quantum 
width of the nucleus, and this contributes to jji2,vv' and higher moments. The other 
contribution is due to transitions between Born-Oppenheimer surfaces. If a nucleus 
starts a trajectory on one particular Born-Oppenheimer energy surface, and there is 
some coupHng to another surface, then that trajectory will split into two trajectories, 
one on each surface. Since the forces experienced on each surface can be different the 
two trajectories can deviate from one another, producing broadening of the total wave 
packet. The first moment (/ti,!^) contains only this contribution, but it appears in all 
other moments as well. 

We now seek to understand the meaning of the moments. First we note that 
Pe, and jl2,uv' are electronic operators (the traces have only been taken over 
the nuclear degrees of freedom). We can therefore take matrix elements with respect 
to electronic states. Suppose we have energy surfaces a characterised by a set of 
electronic states |a). We then have 

Pe,aa =(a|pe|a)= JpaaiR)dR 

Pl^u.aa = (a|Al,iv|") = j Paa{R) (Ru - K) dR 

P2,uv'.aa = {a\il2Mu'\a) = j Paa{R) {Rv ^ Rv){Rv' ~ R„') dR (36) 

where Paa{R) — {aR\p\aR) , and is the ionic density projected onto surface a. Thus 
Pe,aa is just the probability of being on surface a. If we define the mean value of 
some observable Q on surface a to be {a\Q\a) / pe,aa then ^i^v,aa just equals the mean 
position of the nucleus on the surface, measured relative to the mean trajectory R^, 
multiplied by the probability of being on that surface. Similarly, P2,vv',aa gives the 
width of the packet moving on surface a, multiplied by the probability of being on 
that surface, and thus is a measure of the quantum width of the nucleus. 

Thus, in summary, the small amplitude moment expansion corrects the Ehrenfest 
approximation by allowing different trajectories on different energy surfaces, and by 
giving the nuclei a finite quantum width. To make this into a practical method we 
need some efficient way of evaluating the moments fii^v, ti-2,vw' , etc. This is achieved by 
integrating their equations of motion [120]. These equations follow from the Liouville 
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equation (Eq. [T8|l . and the definition of the moments. Thus we have 
dpe 1 



dt 

dAi,t/ 
dt 

dXi,u 
dt 
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(AF^Pe + PeAP^) {^'^^'l^^'^' + Al,.^'^^^') + He{R),\l,u 4(37) 



where = -dHe{R)/dR, 
Trjv|pAA.} and AP^ = 



AF, = - F., i^..' = d^He{R)/dR,dR,', Ai,, = 

Pi/ — Pi/. We can obtain a rather straightforward 

interpretation of these equations by considering matrix elements of the moments with 
respect to the electronic states \a), provided He{R)\a) = ihd\a)/dt. Under these 
conditions from Eq. |S3we obtain 
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dt 
d/ii 
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AF^/ q-qPc Q,Q -|- ^ ^ ^ ( AFjy q,q/ q,/q -t- Pe^CYCt' AFi/(^> "]{38) 

The first equation (for the diagonal elements of the electronic density matrix) tells 
us that transitions between Born-Oppenheimer surfaces are driven by force couplings 
between the surfaces. If we think of p\^v,aa as the position of a nucleus on the surface 
multiplied by the probability of being on the surface, and \\^vMa as the momentum 
of a nucleus on the surface multiplied by the probability of being on the surface (see 
Eq. then the second equation just gives the usual relation between velocity and 
momentum (provided the probability of being on the surface is not changing) . Finally, 
the third equation relates the rate of change of momentum of the nucleus on the surface 
to the force that it is experiencing (the diagonal force term), plus some additional 
corrections derived from the non-adiabatic interactions. Thus we see that the small 
amplitude moment expansion at lower orders is describing classical trajectories on 
multiple Born-Oppenheimer surfaces. 

For a practical implementation of the scheme we need to truncate the infinite 
hierarchy of equations of motion for the moments. If we work with a fixed order of 
expansion for the Hamiltonian (in practice dropping cubic and higher terms in the 
Taylor expansion in Eq. I34II . then for the highest order moments for which we have 
equations of motion there will be additional moments appearing in those equations 
for which there is no corresponding equation of motion. To produce closure we 
therefore need to estimate these higher moments from the lower ones already evaluated. 
For metallic systems, where nuclear wave packets mostly breathe without splitting, 
satisfactory results have been obtained [121,122] by a mean field approximation which 
relates p2,vv' and other second moment operators to products of their own traces and 
Pe. This procedure is a simple case of a general formalism which is currently under 
development. First we reconstruct an approximate density matrix from the known 
moments, then we perform the necessary traces with this density matrix to produce 
the required additional moments. This process is made much easier because of two 



The transfer of energy between electrons and ions in solids 



30 



remarkable theorems. Consider a density matrix for which a Wigner transform has 
been carried out over the nuclear decrees of freedom (see |Appendix B| | , to give pw {RP) 
which is an electronic operator which depends parametrically on the nuclear positions 
and momenta. For one nuclear degree of freedom the first theorem states [216] that 

R''P"'pwiRP)dRdP = (-) ^CrTr|pi?'P'"i?"-'} (39) 

which relates the moments of the Wigner matrix to those of the density matrix. The 
second theorem states that for quadratic Hamiltonians the moments appearing in 
CEID can always be written in the form given by the right-hand side of Eq. EHl In 
short, CEID moments are moments of the Wigner function. This has two immediate 
consequences. First, we do not have to worry about the order of the operators in 
the moments, so that we can completely characterise the moments by the powers of 
position and momentum. Second, since the Wigner function is a function of scalars we 
can approximate this function using standard techniques, rather than attempting the 
much more complex process of approximating functions of operators. For CEID we 
use harmonic oscillator eigenfunctions to expand the Wigner function (an approach 
which has been used in other similar contexts [217]). The final result is that the 
unknown moments can be expressed as a linear combination of the known, with 
constant coefScients§. An alternative approach (that has not been investigated with 
CEID) has been used in a related method [218]. 

To use this method with popular approximate electronic structure methods (such 
as Hartree-Fock or tight binding - and possibly density functional theory, though 
there are additional theoretical problems in this case) it is necessary to reduce the 
N-electron problem into a 1-electron problem (that is, an extension of Hartree-Fock 
theory). This is achieved by taking the A^-electron equations of motion (Eq. ITzjl and 
taking a trace over all the electrons except one. This produces equations of motion 
for 1-electron operators. However, in those equations of motion 2-electron operators 
appear. To reduce the theory to a pure 1-electron theory we need to replace these 2- 
electron operators with suitable functions of 1-electron operators. This is based on the 
density matrix product used in the Hartree-Fock approximation. In matrix notation 
the product is p(12, 1'2') w p{l, l')p{2,2') - p{l,2')p{2, 1'), where the numbers 1 and 2 
refer to the set of indices for electrons 1 and 2 respectively. However, this equation by 
itself is not enough. There are two reasons: first, we need to approximate 2-electron 
moments as well as the electronic density matrix [120]; second, since the electron 
density matrix used in CEID involves a trace over multiple nuclear configurations, we 
need to take an appropriate average of products of single particle density matrices, 
even for the electron density matrix [121]. In both cases the starting point is the same, 
namely to write the moment (q) we are interested in terms of the product of nuclear 
fiuctuations (Q) in the following way 

q = Ttn 
= J Q{R'R)p{RR') AR AR' 

= j Q{R' R)pN{RR')pe[RR') AR AR' (40) 

§ There are also additional correction terms that appear in the equations of motion for the moments 
because of the truncation. These emerge by deriving the equations of motion from an effective 
Lagrangian. This work will be described in a forthcoming paper. 
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where p{RR') = {R\p\R'), Pn{RR') = Tre |/5(i?i?')| and Pe{RR') x pn{RR') ^ 

p{RR'). The operator pe(RR') is an electronic density matrix to which we apply 
the Hartree-Fock approximation. To evaluate the integrals over nuclear coordinates 
we need to make an approximation for the R dependence of pe{RR'). We make a 
Taylor expansion in powers of about Pe{RR), which corresponds to weak dynamic 
coupling between the electrons and nuclei. For details see references [120,121]. Here 
we briefly consider the the result of the above procedure for the electronic density 
matrix: 

p(2)(12; 1'2') « p(i)(ll>W(22') - pi'\l2')pi'H2l') + 

E ^^l' {l^[i'i^M]lm') A^S,(12>K(21')) + . . . (41) 

where pi^^(12; 1'2'), pi^\ll') and ^^^^,(11') are the matrix representations of the two- 
electron density matrix, the one electron density matrix and the one electron first 
moment respectively. The matrix is the inverse of C^J! = {p-2.uiy'}- We 

see that the density matrix is not quite idempotent even for noninteracting electrons, 
which is because there are electron-ion correlations and a dynamical response of the 
electrons to nuclear fiuctuations. This response, described by the second term, screens 
the ion-ion interactions. The resultant dynamical stiffness corrections are essential for 
getting the correct phonon structure and inelastic electron-phonon spectrum [121,122]. 

Historically, CEID was developed to allow electric current induced heating 
in nanoscale devices to be modelled. Therefore, open boundaries have been a 
consideration from the beginning, and have now been implemented in two ways. 
The first starts from Eqs. ESI and EH given above [123,219]. To produce a finite 
system with open boundaries we consider our system as being embedded in an infinite 
environment, and apply our equations of motion to the infinite combined system. 
Clearly we cannot explicitly evolve the electron density matrix and the moments for 
an infinite system, so we make use of the fact that we can write down analytic solutions 
for the evolution of the electron density matrix and the moments provided that the 
Hamiltonian does not vary with time. This allows us to write down closed form 
solutions for the environment which we then couple to the explicit time evolution of the 
system in which we are interested. To produce numerically stable solutions it is found 
necessary to introduce a small amount of damping into the environment. For details 
see references [123,219]. The second scheme, which is currently under development, 
achieves numerical simplifications at the expense of additional approximations. 

So how does this method compare with others? We have already discussed its 
relationship to the Ehrenfest approximation, and clearly it introduces those essential 
features to the theory that allow a proper transfer of energy between electrons and 
nuclei. From Eq. EHl we can see some correspondence with surface hopping (see 
section l3.3.2|l . but with the important difference that different trajectories are able 
to remain coherent with one another as they are treated simultaneously rather than 
independently. A further important feature is that this method reproduces results 
of Fermi's golden rule for the exchange of energy between electrons are nuclear 
provided that the second moment is retained [121]. More broadly, the strengths of 
CEID are that it is not perturbative, it does not invoke the notion of phonons and 
is inherently anharmonic, it avoids classical interpretations of quantum transitions. 
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and it in principle contains the mutual screening of the three interactions (electron- 
electron, electron-nucleus, nucleus-nucleus), while retaining the conceptual framework 
of molecular dynamics. 

CEID is still a rather young method, so a limited range of results has so far 
been produced. It has been applied to Joule heating in atomic wires [120], in which 
case it was found that the first moments {fii.i, and Ai^^) were essential for producing 
the correct heating. The second moments (which involve ARi,ARi,i, AP^AR^r and 
APi^AP^i ) are needed to describe the change in electrical resistance due to increased 
nuclear vibrations [121,122]. Once these were included, it was possible to reproduce 
the spectral signature of the inelastic scattering of electrons by nuclei. 

4. Conclusion and future directions 

Above we have surveyed some of the phenomena resulting from, and the current state 
of theories for understanding, the exchange of energy between electrons and nuclei. 
Accurate modelling of the phenomena is a hard problem because of its intrinsically 
many-body and correlated nature. As a consequence, all the theories are in need of 
further development, as has been indicated in the text. As the different theories have 
different attributes (some are perturbative, others are based on molecular dynamics, 
and so on), they can naturally be applied to different problems. So there is a need for 
development on more than one front. 

The biggest advance is probably going to be in the range of problems to 
be addressed using these techniques. There are now standard phenomena that 
are studied or used to test novel methods, such as heating in nanocontacts, or 
photodissociation, but they are limited in scope. Non-equilibrium phenomena, by 
contrast, are ubiquitous. Certainly the modelling of molecular electronic devices 
and biological molecules associated with non-equilibrium electrons (such as DNA and 
retinal) have already begun to be investigated, and their importance must surely 
increase. Other problems that should benefit from new methods include the evolution 
of radiation damage and inelastic tunnehng spectroscopy carried out with STMs. 
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Appendix A. Electron-phonon Hamiltonians 

There are two common approaches to building Hamiltonians that couple electronic 
and nuclear dynamics (adiabatic states and static lattice states). Even though they 
have been shown to give formally the same answers in lowest order perturbation 
theory [220], there may be practical reasons for choosing one over the other. We 
summarise the basic equations for each case below. Occasionally it is useful to use the 
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diabatic representation of the electronic states. This is obtained from the adiabatic 
representation discussed below by a unitary transform that makes the nuclear kinetic 
energy operator diagonal. 



Appendix A.l. Adiabatic states 

The Hamiltonian for a system of electrons and nuclei we write as H ~ TV + He{R), 
where T/v is the nuclear kinetic energy, He{R) is the Hamiltonian for electrons in 
the field of the nuclei, and R is the position operator for all the nuclei. In the Born- 
Oppenheimer (adiabatic) separation, we assume that we can neglect the nuclear kinetic 
energy to begin with. This then produces the following Schrodinger equation for the 
electrons in the field of the nuclei 

He{R)<^n{R) = S„(i?)$„(i?) (A.l) 

Note that nuclear positions are now represented by numbers i?, and not operators. 
The subscript n indexes the Born-Oppenheimer surfaces. Now that we have a set of 
electronic states <I>„(^), adiabatic nuclear states XnN can be generated by treating 
En{R) as an effective potential for the nuclei, giving 

{Tn + En{R)j XnN = UnNXnN (A. 2) 

where the subscript N indexes the allowed nuclear states, and UnN is the total 
adiabatic energy (electrons and nuclei). Equipped with these states, we can now 
use Fermi's Golden Rule to find transition rates between the product states '^nN = 
XnN^7i, giving the matrix elements MnNn'N' = J di?df 'I'* jyiJ'I'n'W' 



/ dRX*nNXn'N' f dr C&^fjv^n' + ^ ^f^^ f dr 
J J ^ M„Xn'N' J 



(A.3) 



where r is the set of electronic positions, the subscript v indexes the individual nuclear 
coordinates, P^, is the corresponding nuclear momentum operator, and the nuclear 
mass. The energies appearing in Fermi's Golden Rule are UnN- If we assume that 
the energy surfaces on which the nuclei are moving are harmonic, then the nuclear 
wave functions are the simple harmonic oscillator wave functions characterised by a 
frequency, a mass, and an equilibrium position. Further, the energies will satisfy 

UnN = En{Rn,a) +^(nNa + '^j fti^n.a (A.4) 



where Rn,o corresponds to the equihbrium positions of nuclei on surface n, and a is 
an index running over the normal modes which have angular frequencies ujn,a and 
occupancies riNa- The surface dependence of the vibrational frequencies is often 
neglected, that is cj„_q, = uJo,a- 

Appendix A. 2. Static lattice states 

This approach starts with the following partitioning of the Hamiltonian: 

H = fN + H,{R) 

= [Tn + Eo{R) - EoiRo)] + HARq) + [{He{R) - Eo{R)) - (HeiRo) - So(i?o))] 

^ Hn + He{Ro) + HeN (A.5) 
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where Rq is the equihbrium positions of the nuclei in the ground state, and Eq[R) 
is the ground state Born-Oppenheimer energy surface. The electronic states <&„ are 
given by 

7?e(i?o)$„ =-B„$„ (A.6) 

There are two important limits in which we can define the states of the nuclear 
subsystem. When the electron-nuclear coupling HeN is always weak then the reference 
nuclear wavefunctions xn are found from 

HnXn = WnXn (A. 7) 

Fermi's Golden Rule can then be used to find transition rates between the product 
states V^nAT = xn'^u, giving the matrix elements 

MnNn'N' = [En + Wn) 5 N N ' 5nn' + j dRdrX*NKHeNXN'<^n' (A.8) 

In the Hmit of small ionic displacements we can use a linear approximation for the 
interaction Hamiltonian, which is HeN ~ {R — Ro) ■ ^He{Ro), where we have used the 
definition of the fixed sites, namely V£'o(i?o) = 0. Matrix elements of this Hamiltonian 
have the form 

J dRdrX*NKHeNXN'^n' = J dR X*N (R - Ro)XN' ■ J df$;Vi?e(i?0)*n' (A.9) 

The energies appearing in Fermi's Golden Rule are UnN = En + Wn- If the ground 
state Born-Oppenheimer energy surface is harmonic, then the nuclear wave functions 
will be simple harmonic oscillator wave functions centred about the ground state 
equilibrium positions. The nuclear energies will satisfy 

riNa + -z ] fioJa (A. 10) 



where a is an index running over the normal modes which have angular frequencies 
LUa and occupancies riNa- 

The second case corresponds to the electron-nuclear coupHng being weak only 
between energy surfaces, but strong on a surface. In this case the minimum 
energy configurations of different adiabatic energy surfaces are displaced significantly 
from one another. If we treat the electron-nuclear coupling HeN with the 
Hnear approximation we then partition F = ~VHe{Ro) into diagonal (-Fd) 
and off-diagonal (Pod) terms, where / dr (f>* i^£)$„' = 5„„/ J dr $*-?'$„' and 
J df^n^OD^n' = (1 ^ ^riTi') / df^$*j-F$„/. The diagonal Schrodinger equation 
then becomes ^i?Ar + He{Ro) — Fu ■ {R — Ro)j ^nXnN ~ UnN^nXnN, where now the 
nuclear states XnN depend on electronic state. If we make the Harmonic approximation 
for the nuclear Hamiltonian Hn = Tn + ■ K. ■ u, where K is the matrix of spring 
constants and u = R — Rq is the displacement of the nuclei relative to the reference 
positions, then we get 

fN + ^{u~Un)-K-{u~ Un) + En - • K • XnN = UnNXnN (A. 11) 

where the displacement Un is defined by K • ■«„ = J$*FD$„dr. This is just the 
equation of motion for a shifted oscillator whose potential minimum is at u„, and 
whose energy at the minimum is En,o = En — - K-Un. The oscillator wavefunctions 
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are just given by XnN{u) = XnN{u — Un)- The energies appearing in Fermi's Golden 
rule are now 

UnN = E,uQ + ^ (nNa + ^ ) (A. 12) 



2 

and the matrix elements are 

MnNn'N' = ^(1 " Kn') J K^'^n' df • J X*N ~ U„)uXN' {u ^ Un') dR (A.13) 

Appendix B. The Wigner transform 

CEID is formulated in terms of the density matrix, which is then characterised by 
moments of the position and momentum fluctuations of the ions. These moments are 
reminiscent of moments of classical phase space distributions, and so it is natural 
to seek a formal connection. This can be achieved by appealing to the Wigner 
matrix [119,216,221] which is constructed by applying a transformation to the density 
matrix p which is a function of both electronic and nuclear degrees of freedom. Let us 
use a real space representation {\X)) of the TV nuclear degrees of freedom, and leave 
the electronic ones abstract, giving p{X,X') = {X\p\X'). If we make the following 
Hnear combinations, R ^ (X + X')/2 and S ^ X — X' , and carrier out a Fourier 
transformation with respect to S, we get the Wigner matrix 

pw{R, P) = ^ j p{R+ \s. R - \s) cxp(5 • P/ih) dS (B.l) 

The connection between this function and a classical phase space distribution can be 
seen from the following 

(i) If we integrate pw{R,P) over P, we get back the quantum spatial distribution 
function p[R, R). 

(ii) If we integrate pw{R,P) over R, we get the quantum momentum distribution 
function p{P, P). 

(iii) It has an equation of motion similar to that of the classical Liouville equation in 
the classical limit deflned hy h ^ 0, which also corresponds to heavy nuclei [222] 

dpwjR, p) ^ 1 

dt ih 

_ P_ dpwjR, P) 1 j dHejR) dpwjR, P) dpwjR, P) dHejR) 

M dR 2 I ai? dP dR dR 

+ Ojh) 



HejR),pwjR, P) 



Appendix C. Equilibrium and non-equilibrium Green's functions 

In this appendix we give the definition of the different Green's functions, we also 
briefly introduce the technique of the Keldysh time-contour applied to non-equilibrium 
conditions, and derive the corresponding Green's functions and self-energies in the 
presence of interactions. 
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Appendix C.l. General definitions 

Let us consider two operators A and B. The retarded r, advanced a, time-ordered t 
and antichronological time-ordered t Green's function with real time arguments are 
defined as 

G\,B{t,t')^ - m^mA{t),B{t'u), 
G%s{t,t')^m' -t)mtiB{t')]±), 

G^Bit, t') = - '^{Tt {A{t)B{t'))) = -i9{t - t'){A{t)B{t')) ± iOit' - t){B{t')A{t)) , 

G%B{t, t') = - i{Ti {A{t)B{t'))) = -\e{t' - t){A{t)B{t')) ± \6{t - t'){B{t')A{t)) (C.l) 

The average (. . .) is taken over the many-body ground state and A{t) is given in the 
Heisenberg representation. The -I- sign applies when the operators A and B satisfy 
fermion anticommutation relations, and the — sign applies if A and B are boson 
operators. For some applications it is useful to consider Green's functions without 
time ordering. They are the so-called greater > and lesser < Green's functions: 

G%B(i^ t') = -i{A{t)B{t')) and G< s{t, t') = ±\{B{t')A{t)) , 

with the same sign convention as above. All the other Green's functions can be defined 
in terms of these two Green's functions as 

G%B{t, t') = 0{t - t')[G>.,,{t, t') - G%s{t, t')] , 

G%B{t, t') = e{t' - t)[G< s(t, t') - G> s{t, t')] , 

G%Bit,t') = e{t - t')GXs{t,t') + e{t' - t)G%sit,t') . (C.2) 

At (thermal) equilibrium or in a stationary state regime, the Green's functions depend 
on the time difference only {i.e. G\.g{t, t') = G^.g(t—t')), and their Fourier transform 
is dependent on only one energy argument G^{uj). The advanced and retarded Green's 
functions G°''^{uj) contain information about the spectral density of the system, while 
the lesser and greater Green's functions G'^'^{lo) contain information about both the 
spectral density and the occupancy of the system at or out of equilibrium. 

Appendix C.2. Non- equilibrium conditions and time loop contour 

Here we briefly explain the principles of non-equilibrium Green's function [84,85]. 
Consider a many particle system with interactions and/or with a coupling to an 
external driving force (external field). We want to study the system by reducing 
its mathematical description to the calculation of a perturbation series, with the 
hope that later on we can resum some (if not all) the contributions, as is usually 
done in many-body statistical physics [99, 223, 224]. Therefore, we start from 
the non-interacting ground state at the infinitely remote past (where there are no 
interactions) and the interaction V is turned on adiabatically. Then the different 
Green's functions are calculated by going to the interaction picture and evaluating 
the terms of the perturbation expansion series with Wick's theorem. This theorem 
applies to a time-ordered average with respect to a Hamiltonian that is quadratic 
in creation/annihilation operators (for example, a non-interacting Hamiltonian). 
According to Wick's theorem, the average of any product of operators can be found 
by forming all pairs of operators and replacing these by their average. 

For a non-equilibrium (interacting) system, the ground state in the future is not 
known a priori, and we are left with average products which are only partially time 
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ordered. The Keldysh recipe is to introduce a time contour along which the operators 
can be ordered. The time contour Ck contains two branches, the upper (+) and 
the lower (— ) branch. On the upper branch, time starts in the infinitely remote 
past and evolves forwards, then at the turning point (which can be placed at any 
arbitrary time) , one passes onto the lower branch where the system evolves backwards 
in time back to the initially non-interacting starting point aX t = — oo. Then any 

expectation value of products of operator reduces to (0o|T'ca' (^^{t)B(t') . . . Sck^ \4>o) 

where (0o| . . . \4>o) is the average over the non-interacting ground state. The operators 
are given in the interaction picture, i.e. A{t) ~ e'^°*/''A Q-^Hot/h^ ^j^g 
generalization of the time evolution operator (Eq. on the time loop contour 

Sck = ^exp{— i/?i AtV{t)}^ where Tc^ is the time-ordering operator on the 

contour Ck and /^^ At implies integration over Ck- With these definitions, any time 
ordered product can be calculated using the usual rules of many-body perturbation 
theory (Feynman diagrammatic expansion, Wick's theorem, etc.) [99,223,224]. The 
Keldysh recipe is equivalent to reducing the problem to the calculation of averages over 
the non-interacting ground state, which is a great achievement because such averages 
can be calculated exactly in a lot of cases [100]. However, there is a price to pay 
for that: now we have to work with four different Green's functions defined by the 
position of the two times on Ck- When the two times (i, i') are on the same 

branch, the time ordering Tcj^ is equivalent to the standard time ordering: forward 
time ordering on the upper branch and backward time (or anti-time) ordering on the 
lower branch. When {t,t') are on different branches, the time ordering is such that 
any time on the lower branch is always later on the time loop contour Ck than any 
time on the upper branch. 

The electron Green's function defined from the fermion operator \1/ (with 
the definitions in section [Appendix C.l| A ^ ^ and B = *^) G{t,t') = 
—\{Tcii (^(i)^'''^ (<'))) has four Keldysh components on Ck'- for {t^t') on the upper 
branch {+) G{t,t') = -{{Tc^ (*(t+)*t(t'^))) = G\t,t')- for (t,t') on the lower branch 
(-) G{t,t') = -i{Tc^ (*(t_)^'t(i'_))) = G\t,t')- for t on the {+) branch and t' on 
the (-) branch G{t,t') = -\{Tci, (*(t+)^'t(t'_))) = G<{t,t')- and finally for t on the 
(-) branch and t' on the {+) branch G{t,t') = -i{Tc^ {'^{t^)^^t'_^))) G>{t,t'). 
However these Green's functions are not completely independent. They satisfy the 
following relations: G* + G* = G< + G> and G*^ - G° = G> - G< (w ith the retarded 
and advanced Green's functions defined as in section [Appendix C.l| with A = ^ and 
B = ^ft). The self-energy S, associated with the interaction V for example, also has 
four components on the contour Ck- 

It can be shown that the Green's function in the presence of interaction G 
is related to the Green's function in the absence of interaction Go via the usual 
Dyson equation G{t,t') = Go{t,t') + J^^, dtidt2 Go(t, ti)S(ii, t2)G(t2, i') where the 
time integrals are taken over Ck- The Dyson equation can then be re-expressed 
using only integrals over the real-time axis by introducing the different components of 
G and E on the contour Ck- By using the relation between the different Green's 
functions, we finally obtain a Dyson-like equation for the advanced and retarded 
(non-equihbrium) Green's functions: G'''° = Gg " -I- Gq ''S'''"G'''" and another kinetic 
equation for the (non-equihbrium) lesser and greater Green's functions: G^'^ = 
(1 + G''E'^)G;5^'>(1 + S''G'') + G'-S<'>G'^. In these equations, the products GS or SG 
imply time integration over the real axis, i.e. dt. Similar results can be derived 
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for the phonon Green's functions {i.e. when A and B are phonon field operators). 

In principle, we can now calculate exactly the non-equilibrium properties of a 
many-body interacting system by determining self-consistently the different Green's 
functions and self-energies (knowing that the self-energies are functionals of the 
different Green's functions themselves). 



Appendix C.3. Electron and phonon Green's functions 

Appendix C.3.1. Non-interacting systems For quadratic Hamiltonians {i.e. non- 
interacting particles), we can calculate exactly the different Green's functions. Let 
us start with electrons and consider the electronic Hamiltonian Hq = X]n^nc!j^c„. 
The operator cj^ (c„) creates (annihilates) an electron in the rt-th electronic state 
with energy e„; n labels either the eigenstates of a finite size system or the fc-states 
of a (l,2,3)-dimensional periodic system. With the definitions given in the previous 
sections and taking A = c„ and B = A\ we find the following Green's functions in 
the energy representation: 

G"'"(w) = with ?/ ^ 0+ , 

G<{ij) = 27ri(4c„) 6{l, - e„) , G>{u;) = -27ri(l - c„)) 6{iJ - e„) . (C.3) 

The average (cJjC„) gives the occupation number of the state n, and in the 
thermodynamic Hmit it is given by the Fermi-Dirac distribution /(£«) = i/(e'^(^"~'^) + 
1) where ^ is chemical potential and f3 = l/fc^T. 

For phonons with a quadratic Hamiltonian Hq = J2\ f'^x{o-\c-\ + 1/2) where a\ 
{a\) creates (annihilates) a quantum of energy hojx, we find the following Green's 
functions (denoted by the letter D) defined from A = a\ + ax and B = A^: 

i?S'"(w) = : ^with?7^0+, 

D<^>{lu)= -27ri{{nx)S{ujTLOx) + {l + {nx)) S{lu±lox)}, (C.4) 

where nx ~ a\ax. In the thermodynamic Hmit, {nx) is given by the Bose-Einstein 
distribution N{u;) = l/(e'^'^ - 1). 



Appendix C.3. 2. Coupling to reservoirs If the (finite size) system of interest is 
connected to M other subsystems (generally much larger) which act as either thermal 
reservoirs or particle reservoirs, then we can calculate the full Green's function of 
the connected primary system (or central region) from the Green's function of the 
isolated parts of the system. This is done by solving the Dyson equation and the 
quantum kinetic equation for the non-equilibrium Green's functions. As an example, 
we consider that the coupling to the M subsystems is given by the Hamiltonian 
matrices Va corresponding to electron hopping between the primary system and the 
M other subsystems. When neglecting the interactions between particles, the self- 
energies Y!i^=(^^°-^<^>) entering the Dyson equations for the electron Green's functions 

of the central region (coupled to the M reservoirs) are = X^qLi where 
^a('^) = ^agaai^Wa and g^^ is the corresponding Green's function of the isolated 
a-th subsystem fEn. lC.3|l . 

For a central region consisting only of phonon modes A coupled to another phonon 
bath (set of phonon of frequency ujf}) via some coupling contants Up, one can derive 
the full phonon Green's function of the central region as 
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where fg" is the bare phonon Green's function given in En. IC.4I and the self-energy 
n^'° arises from the coupHng of the modes A to the modes f3. In the simplest case, 
n^'° can be approximated by n^'"(a;) oc iJ2i3 ~ ^fi) [225]. 

In more realistic systems, the central region needs to be described by a 
Hamiltonian that also includes interaction between electrons or between electrons and 
phonons. In the next section, the case of electron-phonon interaction is considered in 
detail. 



Appendix C.3.3. Self-energies for electron-phonon interaction A general form for 
the linear electron-phonon (e-ph) interaction Hamiltonian is as follows: -f^gpjj = 

n m l>'nm{a\ + a\)c\^Cm where 7Anm are the matrix elements of the e-ph coupling 
matrix 7_)^. For such an interaction, one needs to include another contribution S^p/j 

in the electron self-energies = X^aLi^S- Similarily, one has to include the 
contribution due to the e-ph interaction in the phonon self-energies 11^. 

As mentioned in section 13.2.31 within the self-consistent Born approximation, 
the lowest-order perturbation expansion for the interaction is used to determine the 
self-energies and then the non-interacting Green's functions is substituted by the full 
Green's functions of the system. As an example, the contribution from the Fock- 
diagrams to the electron self-energy is: {t,t') oc i X^a -^(^' ^')7a^(^j ^')7a ^-^d 
t' being on the contour Ck)- Using Langreth's rules for products of operators on 
Ck [226], one gets the different contributions to the self-energy in the following 
form: Yj''J^^[uj) oc i j dijj' D'-^[uj — uj')'^)^Gy{uj')'^x^ where x,y = (r,a,<,>). For 

F r 

example, the retarded Fock part of the electron self-energy is given by: S = 
^eph ~'~ ^eph ^eph' exact expressions for the different self-energies can be 

found in Refs. [103,104,115,116,227]. The contribution to the phonon self-energies 
can be calculated to second order in the e-ph coupling by considering diagrams 
for electron-hole excitations. The polarisation IIa can be obtained from \lx{t,t') oc 
\Tx{^xG{t,t')'^x^{t' ^t)} where the trace runs over the electron states n. Once again, 
one has to use Langreth's rules to obtain the different contributions li'^^'^'^ {to) which 
can be found in Refs. [106, 107, 109, 115]. And finally, we solve the problem by 
calculating the different electron and phonon full Green's functions G°''''^'''(u;) and 
£)a,r,<,>^^^ from the coupled Dyson and quantum kinetic equations, with the self- 
energies being functional in the following form: Yf^'^^'^ {uj) = S[{G'^(a;)}, 

and IiY^<^''{^) = n[{G-(c.)}, {D-{uo)}]. 
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